##This is the code for figures and statistics from Barnes et al. (2024): phosphorus composition along a burn severity gradient ##Contact: Morgan E Barnes ##Make your environemnt

Set Up

Load libraries

require(pacman)
## Loading required package: pacman
pacman::p_load(tidyverse,
              dplyr,plyr,ggplot2,
              lubridate, tidyr, tidyverse,
              scales, cowplot, gridExtra, RColorBrewer,
              wesanderson,readxl,openxlsx, stringr, hablar,
              janitor, stringr, naniar, forstringr, reshape2, 
              emmeans, rstatix, ggpubr, Hmisc, corrplot, ggrepel, 
              stringr, lme4, lmerTest, car, performance, vegan, psych, lavaan)
Confirm working directory:
## [1] "/Users/barn772/OneDrive - PNNL/Documents/GitHub/rcsfa-RC3-BSLE_P/scripts"

##Load data The datasets and metadata needed for P BSLE Gradient manuscript are: 1) Chemistry datasets solids and liquids w/ all metadata that can be used directly from data package: This includes: a) Total elemental composition of solids b) Total elemental composition of leachate particulate and filtered fractions c) Molybdate reactive P of leachate filtered fraction d) XANES of solids e) NMR of solids

  1. Chemistry data that needs post-data package processing after download from BSLE Data Package v3; this was done and provided within the BSLE P Manuscript Data Package:
  1. XANES of solids (spectra of samples and reference compounds and linear combination fits)
  2. NMR of solids (deconvolution)

Download the BSLE Data Package

Until we can get the API to work, do this manually and store on your local. I recommend storing inside your github folder and gitignoring the subfolder with the data in it by adding this to your .gitignore: #data folder data/

BSLE v3 data package URL: https://data.ess-dive.lbl.gov/view/doi:10.15485/1894135 BSLE P manuscript data package URL: INSERT HERE!!! fix

Bring in the data package data:

#metadata
Meta <- read_csv("../data/BSLE_Data_Package_v3/v1_BSLE_Metadata_and_Protocols/BSLE_Burn_and_Laboratory_Metadata.csv") %>%
  dplyr::rename(Sample_ID = Sample_Name,
                Leachate_vol_mL = Leachate_Volume,
                Char_leached_g = Char_Leached) %>%
  naniar::replace_with_na_all(condition = ~.x ==-9999) %>%
  naniar::replace_with_na_all(condition = ~.x == "N/A") %>%
  mutate(Burn_Severity=replace(Burn_Severity, Burn_Severity=='unburned', 'Unburned'))
## Rows: 234 Columns: 26
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (19): Burn_Code, Parent_ID, Sample_Name, Feedstock_Species, Land_Coverag...
## dbl  (7): Max_Burn_Temp, Dry_Fuel_Before_Burn_Weight, Approx_Char_After_Burn...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Meta<-Meta[ c(2,3,5,8,21,25,26,6,20,22,23)]
Meta$Burn_Severity=factor(Meta$Burn_Severity, levels= c("Unburned", "Low", "Moderate", "High"))

#chemical data from BSLE v3 data package
ICP_Solid <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_ICP_Solid.csv")
## New names:
## Rows: 33 Columns: 12
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (12): #Columns, 11, ...3, ...4, ...5, ...6, ...7, ...8, ...9, ...10, ......
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
ICP_Leachate <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_ICP.csv")
## New names:
## Rows: 128 Columns: 12
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (12): #Columns, 11, ...3, ...4, ...5, ...6, ...7, ...8, ...9, ...10, ......
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
MBD <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_MBD.csv")
## New names:
## Rows: 71 Columns: 5
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (5): #Columns, 4, ...3, ...4, ...5
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
ICP_PNMR <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_ICP_PNMR.csv")
## New names:
## Rows: 33 Columns: 20
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (20): #Columns, 19, ...3, ...4, ...5, ...6, ...7, ...8, ...9, ...10, ......
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
npoctdn <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/v2_BSLE_NPOC_TN.csv")
## New names:
## Rows: 328 Columns: 6
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (6): #Columns, 5, ...3, ...4, ...5, ...6
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
SolidCN <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_CN.csv")
## New names:
## Rows: 83 Columns: 6
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (6): #Columns, 5, ...3, ...4, ...5, ...6
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
pH <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_pH.csv")
## New names:
## Rows: 247 Columns: 5
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (5): #Columns, 4, ...3, ...4, ...5
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
NMR_spectra <- read_csv("../data/BSLE_Data_Package_v3/v3_BSLE_Data/BSLE_31P-NMR.csv")
## Rows: 374667 Columns: 58
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (58): Chemical_Shift_parts_per_million, BSLE_0001-solid, BSLE_0002-solid...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#chemical data from BSLE P manuscript data package
NMR <- read.csv("../data/Barnes_2023_BSLE_P_Gradient_Manuscript_Data_Package/NMR_Deconvolutions.csv") %>%
  filter(Degradation_Corrected == 'Yes')
XANES <- read.csv("../data/Barnes_2023_BSLE_P_Gradient_Manuscript_Data_Package/P_XANES_Samples/P-XANES_LCF.csv")  
XANES_spectra <- read_csv("../data/Barnes_2023_BSLE_P_Gradient_Manuscript_Data_Package/P_XANES_Samples/xanes_example_spectra.csv") #this file was generated using XANES_Spectra_Merging script 
## New names:
## Rows: 326 Columns: 9
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," dbl
## (9): ...1, eV, BSLE_0013-solid_P-XANES.xmu.nor, BSLE_0007-solid_P-XANES....
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`

Clean up data functions:

essdive_template_cleaning <- function(dataframe){
  cleaned <- dataframe %>% 
  janitor::row_to_names(row_number = which(dataframe[,1] == "Field_Name")) %>% #this assumes that this name is somewhere in column #1. That follows some (most?) ESS-DIVE templates.
 slice(which(Field_Name == "#Start_Data"):which(Field_Name == "#End_Data")) %>%
 select(!Field_Name)
  }

noreporting_format_cleaning <- function(dataframe){
  cleaned <- dataframe %>% 
  janitor::row_to_names(row_number = which(dataframe[,1] == "Date")) #this assumes that this name is somewhere in column #1. This works for the thermocouple datasets
  }
pH_clean <- essdive_template_cleaning(pH)  %>%
  dplyr::rename(Sample_ID = Sample_Name,
                pH = '00403_pH') %>%
   mutate_at(c(3), as.numeric) %>%
  separate(col=Sample_ID, into = c("Sample_ID", "FilterSize"), sep="-")  %>%
  filter(grepl('unfilt', FilterSize)) %>%
  select(Sample_ID, pH)

#wrote modifies so -9999 is NA, below LOD is NA, and if above or below standard curve then I pull out those values
npoctdn_clean <- essdive_template_cleaning(npoctdn) %>%
  dplyr::rename(NPOC_mg_C_per_L = `00681_NPOC_mg_per_L_as_C`,
         TN_mg_N_per_L = `00602_TN_mg_per_L_as_N`,
         Sample_ID = Sample_Name) %>%
  select(Sample_ID, NPOC_mg_C_per_L, TN_mg_N_per_L) %>%
  filter(grepl('filt0.7', Sample_ID)) %>%
  mutate(Sample_ID = stringr::str_remove(Sample_ID ,"-filt0.2"),
         NPOC_mg_C_per_L = case_when(NPOC_mg_C_per_L == -9999 ~ NA,
                              str_detect(NPOC_mg_C_per_L, "LOD") ~ NA,
                              str_detect(NPOC_mg_C_per_L, "Below_0.1_ppm") ~
                                      str_extract_part(NPOC_mg_C_per_L, "_ppm_Below_0.1_ppm", before = TRUE),
                               TRUE ~ NPOC_mg_C_per_L),
         TN_mg_N_per_L = case_when(TN_mg_N_per_L == -9999 ~ NA,
                              str_detect(TN_mg_N_per_L, "LOD") ~ NA,
                              str_detect(TN_mg_N_per_L, "Below_0.1_ppm") ~ 
                                        str_extract_part(TN_mg_N_per_L, "_ppm_Below_0.1_ppm", before = TRUE),
                               TRUE ~ TN_mg_N_per_L)) %>%
  mutate(NPOC_mg_C_per_L = case_when(str_detect(NPOC_mg_C_per_L, "NPOC_") ~ 
                                     str_extract(NPOC_mg_C_per_L, "(\\d+\\.\\d+)"),
                                     NPOC_mg_C_per_L < 0 ~ NA,
                                     TRUE ~ NPOC_mg_C_per_L),
         TN_mg_N_per_L = case_when(str_detect(TN_mg_N_per_L, "TN_") ~ 
                                     str_extract(TN_mg_N_per_L, "(\\d+\\.\\d+)"),
                                     TN_mg_N_per_L < 0 ~ NA, 
                                     TRUE ~ TN_mg_N_per_L)) %>%
  mutate(NPOC_mg_C_per_L = as.numeric(NPOC_mg_C_per_L),
         TN_mg_N_per_L = as.numeric(TN_mg_N_per_L)) %>%
  separate(col=Sample_ID, into = c("Sample_ID", "FilterSize"), sep="-") 

SolidCN_clean <- essdive_template_cleaning(SolidCN) %>%
  dplyr::rename(C_weight_per = `01463_C_percent`,
         N_weight_per = `01472_N_percent`,
         Parent_ID = Sample_Name) %>%
  select(Parent_ID, C_weight_per,  N_weight_per) %>%
  filter(grepl('solid', Parent_ID)) %>%
  mutate(Parent_ID = stringr::str_remove(Parent_ID ,"-solid"),
         C_weight_per = case_when(C_weight_per == -9999 ~ NA,
                                  str_detect(C_weight_per, "LOD") ~ NA,
                                     TRUE ~ C_weight_per),
         N_weight_per = case_when(N_weight_per == -9999 ~ NA,
                                   str_detect(N_weight_per, "LOD") ~ NA,
                                   TRUE ~ N_weight_per)) %>%
   mutate(C_weight_per = as.numeric(C_weight_per),
         N_weight_per = as.numeric(N_weight_per))

#negative values were replaced with NA, below LOD were replaced with NA, and above/below standard curve were replaced with the value
ICP_Solid_clean<- essdive_template_cleaning(ICP_Solid) %>%
  dplyr::rename(Solid_Ca_mg_kg = '52718_Total_Ca_mg_per_kg',
                Solid_Mg_mg_kg = '52719_Total_Mg_mg_per_kg',
                Solid_P_mg_kg = 'Total_P_mg_per_kg',
                Solid_K_mg_kg = '52720_Total_K_mg_per_kg',
                Solid_Na_mg_kg = '52721_Total_Na_mg_per_kg',
                Solid_S_mg_kg = '01466_Total_S_mg_per_kg',
                Solid_Al_mg_kg = '01462_Total_Al_mg_per_kg',
                Solid_Fe_mg_kg = '01464_Total_Fe_mg_per_kg',
                Parent_ID = 'Sample_Name') %>%
  mutate(Solid_Na_mg_kg = case_when(
                Solid_Na_mg_kg == -9999 ~ NA,
                str_detect(Solid_Na_mg_kg, "Negative_Value") ~ NA, 
                str_detect(Solid_Na_mg_kg, "Below_213.43_LOD") ~ NA, 
                str_detect(Solid_Na_mg_kg, "TotNa_Below") ~ 
                str_extract(Solid_Na_mg_kg, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
                str_detect(Solid_Na_mg_kg, "TotNa_Above") ~ 
                str_extract(Solid_Na_mg_kg, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
                TRUE ~ Solid_Na_mg_kg)) %>%
   mutate_at(c(4:11), as.numeric) %>%
  separate(col=Parent_ID, into = c("Parent_ID", "FilterSize"), sep="-")  %>%
  filter(grepl('solid', FilterSize)) %>%
  select(Parent_ID, Solid_Ca_mg_kg,  Solid_Mg_mg_kg, Solid_P_mg_kg, Solid_K_mg_kg, Solid_Na_mg_kg, Solid_S_mg_kg, Solid_Al_mg_kg, Solid_Fe_mg_kg)
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Methods_Deviation =
##   .Primitive("as.double")(Methods_Deviation)`.
## Caused by warning:
## ! NAs introduced by coercion
#negative values were replaced with NA, below LOD were replaced with NA, and above/below standard curve were replaced with the value
ICP_Leachate_clean <- essdive_template_cleaning(ICP_Leachate) %>%
  dplyr::rename(Leachate_Ca_mg_L = 'Total_Ca_mg_per_L',
                Leachate_Mg_mg_L = 'Total_Mg_mg_per_L',
                Leachate_P_mg_L = 'Total_P_mg_per_L',
                Leachate_K_mg_L = 'Total_K_mg_per_L',
                Leachate_Na_mg_L = 'Total_Na_mg_per_L',
                Leachate_S_mg_L = 'Total_S_mg_per_L',
                Leachate_Al_mg_L = 'Total_Al_mg_per_L',
                Leachate_Fe_mg_L = 'Total_Fe_mg_per_L',
                Sample_ID = 'Sample_Name') %>%
  mutate(Leachate_Na_mg_L = case_when(
                Leachate_Na_mg_L == -9999 ~ NA,
                str_detect(Leachate_Na_mg_L, "Negative") ~ NA, 
                str_detect(Leachate_Na_mg_L, "LOD") ~ NA, 
                str_detect(Leachate_Na_mg_L, "TotNa_Below") ~ 
                str_extract(Leachate_Na_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
                str_detect(Leachate_Na_mg_L, "TotNa_Above") ~ 
                str_extract(Leachate_Na_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
                TRUE ~ Leachate_Na_mg_L)) %>%
  mutate(Leachate_Al_mg_L = case_when(
                Leachate_Na_mg_L == -9999 ~ NA,
                str_detect(Leachate_Al_mg_L, "Negative") ~ NA, 
                str_detect(Leachate_Al_mg_L, "LOD") ~ NA, 
                str_detect(Leachate_Al_mg_L, "TotAl_Below") ~ 
                str_extract(Leachate_Al_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
                str_detect(Leachate_Al_mg_L, "TotAl_Above") ~ 
                str_extract(Leachate_Al_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
                TRUE ~ Leachate_Al_mg_L)) %>%
  mutate(Leachate_Fe_mg_L = case_when(
                Leachate_Fe_mg_L == -9999 ~ NA,
                str_detect(Leachate_Fe_mg_L, "Negative") ~ NA, 
                str_detect(Leachate_Fe_mg_L, "LOD") ~ NA, 
                str_detect(Leachate_Fe_mg_L, "TotAl_Below") ~ 
                str_extract(Leachate_Fe_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
                str_detect(Leachate_Fe_mg_L, "TotAl_Above") ~ 
                str_extract(Leachate_Fe_mg_L, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
                TRUE ~ Leachate_Fe_mg_L))  %>%
  mutate_at(c(3:10), as.numeric) %>%
  select(Sample_ID, Leachate_Ca_mg_L,  Leachate_Mg_mg_L, Leachate_P_mg_L, Leachate_K_mg_L, Leachate_Na_mg_L, Leachate_S_mg_L, Leachate_Al_mg_L, Leachate_Fe_mg_L) %>%
  filter(grepl('filt', Sample_ID)) %>%
  separate(col=Sample_ID, into = c("Sample_ID", "FilterSize"), sep="-") 

ICP_Leachate_Filt_clean<- ICP_Leachate_clean %>%
                filter(str_detect(FilterSize, "filt0.")) %>%
  dplyr::rename(Leachate_Ca_mg_L_filt0.7 = 'Leachate_Ca_mg_L',
                Leachate_Mg_mg_L_filt0.7 = 'Leachate_Mg_mg_L',
                Leachate_P_mg_L_filt0.7 = 'Leachate_P_mg_L',
                Leachate_K_mg_L_filt0.7 = 'Leachate_K_mg_L',
                Leachate_Na_mg_L_filt0.7 = 'Leachate_Na_mg_L',
                Leachate_S_mg_L_filt0.7 = 'Leachate_S_mg_L',
                Leachate_Al_mg_L_filt0.7 = 'Leachate_Al_mg_L',
                Leachate_Fe_mg_L_filt0.7 = 'Leachate_Fe_mg_L') %>%
  select(Sample_ID, Leachate_Ca_mg_L_filt0.7,  Leachate_Mg_mg_L_filt0.7, Leachate_P_mg_L_filt0.7, Leachate_K_mg_L_filt0.7, Leachate_Na_mg_L_filt0.7, Leachate_S_mg_L_filt0.7, Leachate_Al_mg_L_filt0.7, Leachate_Fe_mg_L_filt0.7)

ICP_Leachate_Unfilt_clean<- ICP_Leachate_clean %>%
                filter(str_detect(FilterSize, "unfilt")) %>%
  dplyr::rename(Leachate_Ca_mg_L_unfilt = 'Leachate_Ca_mg_L',
                Leachate_Mg_mg_L_unfilt = 'Leachate_Mg_mg_L',
                Leachate_P_mg_L_unfilt = 'Leachate_P_mg_L',
                Leachate_K_mg_L_unfilt = 'Leachate_K_mg_L',
                Leachate_Na_mg_L_unfilt = 'Leachate_Na_mg_L',
                Leachate_S_mg_L_unfilt = 'Leachate_S_mg_L',
                Leachate_Al_mg_L_unfilt = 'Leachate_Al_mg_L',
                Leachate_Fe_mg_L_unfilt = 'Leachate_Fe_mg_L') %>%
  select(Sample_ID, Leachate_Ca_mg_L_unfilt,  Leachate_Mg_mg_L_unfilt, Leachate_P_mg_L_unfilt, Leachate_K_mg_L_unfilt, Leachate_Na_mg_L_unfilt, Leachate_S_mg_L_unfilt, Leachate_Al_mg_L_unfilt, Leachate_Fe_mg_L_unfilt)


#negative values were replaced with zero, below LOD were replaced with NA, and above/below standard curve were replaced with the value (value is corrected for dilution and subtracted by unreactive P to correct for color interference)
#only used concentrations above the high standard if within 15% of it
MBD_Leachate_clean <- essdive_template_cleaning(MBD) %>%
  dplyr::rename(Sample_ID = 'Sample_Name',
                MBD_P_mg_per_L_filt0.7 = 'MBD_P_mg_per_L') %>%
  mutate(MBD_P_mg_per_L_filt0.7 = case_when(
                MBD_P_mg_per_L_filt0.7 == -9999 ~ NA,
                str_detect(MBD_P_mg_per_L_filt0.7, "negative") ~ NA, 
                str_detect(MBD_P_mg_per_L_filt0.7, "LOD") ~ NA, 
                str_detect(Sample_ID, "50A") ~ NA, #hihger than 15% of high standard
                str_detect(MBD_P_mg_per_L_filt0.7, "MBD_Below") ~ 
                str_extract(MBD_P_mg_per_L_filt0.7, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
                str_detect(MBD_P_mg_per_L_filt0.7, "MBD_Above") ~ 
                str_extract(MBD_P_mg_per_L_filt0.7, "\\d+\\.?\\d*(?=_ppm_Final_Corrected)"),
                TRUE ~ MBD_P_mg_per_L_filt0.7)) %>%
  mutate_at(c(2), as.numeric) %>%
  filter(grepl('filt', Sample_ID)) %>%
  separate(col=Sample_ID, into = c("Sample_ID", "FilterSize"), sep="-") %>%
  select(Sample_ID, MBD_P_mg_per_L_filt0.7) %>%
  mutate_at(c(2), as.numeric) %>%
  filter(Sample_ID != "BSLE_0073C") #remove 73C because below limit of quantification
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Material = .Primitive("as.double")(Material)`.
## Caused by warning:
## ! NAs introduced by coercion
ICP_PNMR_clean <- essdive_template_cleaning(ICP_PNMR) %>%
  dplyr::rename(Parent_ID = 'Sample_Name') %>%
  mutate_at(c(2:8), as.numeric) %>%
  filter(grepl('solid', Parent_ID)) %>%
  separate(col=Parent_ID, into = c("Parent_ID", "FilterSize"), sep="-") %>%
  select(Parent_ID, Total_Extractable_Ca_mg_per_kg,  Total_Extractable_Mg_mg_per_kg, Total_Extractable_P_mg_per_kg, Total_Extractable_K_mg_per_kg,  Total_Extractable_S_mg_per_kg, Total_Extractable_Al_mg_per_kg, Total_Extractable_Fe_mg_per_kg) %>%
  mutate_at(c(2:8), as.numeric) 
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Material = .Primitive("as.double")(Material)`.
## Caused by warning:
## ! NAs introduced by coercion
XANES_clean <- XANES %>%
  select(-where(~all(is.na(.)))) %>%
  mutate(XANES_Po = rowSums(select(., contains("Po")), na.rm = TRUE)) %>%
  mutate(XANES_Pi = rowSums(select(., contains("Pi")), na.rm = TRUE)) %>%
  mutate(XANES_Fe = rowSums(select(., contains("Fe")), na.rm = TRUE)) %>%
  mutate(XANES_Na = rowSums(select(., contains("Pi_Na") | contains("Po_Na")), na.rm = TRUE)) %>%
  mutate(XANES_Ca = rowSums(select(., contains("Ca") | contains("Apatite")), na.rm = TRUE)) %>%
  mutate(XANES_Al = rowSums(select(., contains("Al")), na.rm = TRUE)) %>%
  mutate(XANES_Mg = rowSums(select(., contains("Mg") | contains("NH4_Mg")), na.rm = TRUE)) %>%
  mutate(XANES_K = rowSums(select(., contains("K")), na.rm = TRUE)) %>%
  mutate(XANES_Pi_Fe = rowSums(select(., contains("Pi_Fe")), na.rm = TRUE)) %>%
  mutate(XANES_Pi_Na = rowSums(select(., contains("Pi_Na")), na.rm = TRUE)) %>%
  mutate(XANES_Pi_Ca = rowSums(select(., contains("Pi_Ca") | contains("Apatite")), na.rm = TRUE)) %>%
  mutate(XANES_Pi_Al = rowSums(select(., contains("Pi_Al")), na.rm = TRUE)) %>%
  mutate(XANES_Pi_Mg = rowSums(select(., contains("Pi_Mg") | contains("NH4_Mg")), na.rm = TRUE)) %>%
  mutate(XANES_Pi_K = rowSums(select(., contains("Pi_K")), na.rm = TRUE)) %>%
  mutate(XANES_Po_Fe = rowSums(select(., contains("Po_Fe")), na.rm = TRUE)) %>%
  mutate(XANES_Po_Na = rowSums(select(., contains("Po_Na")), na.rm = TRUE)) %>%
  mutate(XANES_Po_Ca = rowSums(select(., contains("Po_Ca")), na.rm = TRUE)) %>%
  mutate(XANES_Po_Al = rowSums(select(., contains("Po_Al")), na.rm = TRUE)) %>%
  mutate(XANES_Po_Mg = rowSums(select(., contains("Po_Mg")), na.rm = TRUE)) %>%
  mutate(XANES_Po_K = rowSums(select(., contains("Po_K")), na.rm = TRUE)) %>%
  separate(col=Sample_Name, into = c("Parent_ID", "FilterSize"), sep="-") %>%
  filter(str_detect(FilterSize, "solid"))

NMR_clean <- NMR %>%
  separate(col=Sample_Name, into = c("Parent_ID", "FilterSize"), sep="-") 
data <- Meta %>%
  full_join(SolidCN_clean, by = 'Parent_ID') %>%
  full_join(ICP_Solid_clean, by = 'Parent_ID') %>%
  full_join(ICP_PNMR_clean, by = 'Parent_ID') %>%
  full_join(npoctdn_clean, by = 'Sample_ID') %>%
  full_join(ICP_Leachate_Filt_clean, by = 'Sample_ID') %>%
  full_join(ICP_Leachate_Unfilt_clean, by = 'Sample_ID') %>%
  full_join(MBD_Leachate_clean, by = 'Sample_ID') %>%
  full_join(pH_clean, by = 'Sample_ID') %>%
  mutate_at(c(10:48), as.numeric) %>%
  filter(!is.na(Solid_P_mg_kg))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `FilterSize = .Primitive("as.double")(FilterSize)`.
## Caused by warning:
## ! NAs introduced by coercion
data$Burn_Severity=factor(data$Burn_Severity, levels= c("Unburned", "Low", "Moderate", "High"))

#create a dataframe for only solids figures and stats to prevent having reps count as extra samples; add NMR and XANES data
data_solids <- data %>%
  filter(str_detect(Sample_ID, "A")) %>%
  select(1:28) %>%
  left_join(NMR_clean) %>%
  left_join(XANES_clean)
## Joining with `by = join_by(Parent_ID)`
## Joining with `by = join_by(Parent_ID, FilterSize)`
data_solids$Burn_Severity=factor(data_solids$Burn_Severity, levels= c("Unburned", "Low", "Moderate", "High"))

#keep XANES separate too
XANES_clean<- Meta %>%
  full_join(XANES_clean, by = 'Parent_ID') %>%
  filter(str_detect(Sample_ID, "A")) %>%
  mutate(Burn_Treatment = as.character(Burn_Treatment)) %>%
  filter(Burn_Treatment %in% c("burn table" , "raw")) %>%
  filter(Land_Coverage_Category %in% c("Douglas-fir forest" , "Sagebrush shrubland")) %>%
  filter(!is.na(FilterSize))

####Additional Calcs

#NMR extraction efficiency
data_solids<-data_solids %>%
  mutate(NMR_EE_P_perc=Total_Extractable_P_mg_per_kg/Solid_P_mg_kg*100)

#convert C and N of solids from percent to mg/kg
data_solids<-data_solids %>%
  mutate(Solid_C_mg_kg = (C_weight_per*10000),
  Solid_N_mg_kg = (N_weight_per*10000))


#convert to P relative to solid mass and water volume for Molybdate
data<-data %>%
  mutate(MBD_P_mg_kgchar = MBD_P_mg_per_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) 

#normalize by P concentration in solid char for Molybdate
data<-data %>%
  mutate(MBD_P_percentofsolid = MBD_P_mg_kgchar/Solid_P_mg_kg*100)

#calculate particulate concentrations as the difference between filtered and unfiltered
data<-data %>%
  mutate(Leachate_Ca_mg_L_part = Leachate_Ca_mg_L_unfilt - Leachate_Ca_mg_L_filt0.7) %>%
  mutate(Leachate_Mg_mg_L_part = Leachate_Mg_mg_L_unfilt - Leachate_Mg_mg_L_filt0.7) %>%
  mutate(Leachate_P_mg_L_part = Leachate_P_mg_L_unfilt - Leachate_P_mg_L_filt0.7) %>%
  mutate(Leachate_K_mg_L_part = Leachate_K_mg_L_unfilt - Leachate_K_mg_L_filt0.7) %>%
  mutate(Leachate_Na_mg_L_part = Leachate_Na_mg_L_unfilt - Leachate_Na_mg_L_filt0.7) %>%
  mutate(Leachate_S_mg_L_part = Leachate_S_mg_L_unfilt - Leachate_S_mg_L_filt0.7) %>%
  mutate(Leachate_Al_mg_L_part = Leachate_Al_mg_L_unfilt - Leachate_Al_mg_L_filt0.7) %>%
  mutate(Leachate_Fe_mg_L_part = Leachate_Fe_mg_L_unfilt - Leachate_Fe_mg_L_filt0.7)

#convert to P relative to solid mass and water volume for the soluble and particulate fractions
data<-data %>%
  mutate(Leachate_Ca_mg_kgchar_part = Leachate_Ca_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Mg_mg_kgchar_part  = Leachate_Mg_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_P_mg_kgchar_part  = Leachate_P_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_K_mg_kgchar_part  = Leachate_K_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Na_mg_kgchar_part  = Leachate_Na_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_S_mg_kgchar_part  = Leachate_S_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Al_mg_kgchar_part  = Leachate_Al_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Fe_mg_kgchar_part  = Leachate_Fe_mg_L_part*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Ca_mg_kgchar_filt0.7 = Leachate_Ca_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Mg_mg_kgchar_filt0.7 = Leachate_Mg_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_P_mg_kgchar_filt0.7 = Leachate_P_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_K_mg_kgchar_filt0.7 = Leachate_K_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Na_mg_kgchar_filt0.7 = Leachate_Na_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_S_mg_kgchar_filt0.7 = Leachate_S_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Al_mg_kgchar_filt0.7 = Leachate_Al_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Fe_mg_kgchar_filt0.7 = Leachate_Fe_mg_L_filt0.7*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Ca_mg_kgchar_unfilt = Leachate_Ca_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Mg_mg_kgchar_unfilt  = Leachate_Mg_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_P_mg_kgchar_unfilt = Leachate_P_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_K_mg_kgchar_unfilt  = Leachate_K_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Na_mg_kgchar_unfilt  = Leachate_Na_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_S_mg_kgchar_unfilt  = Leachate_S_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Al_mg_kgchar_unfilt  = Leachate_Al_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000) %>%
  mutate(Leachate_Fe_mg_kgchar_unfilt = Leachate_Fe_mg_L_unfilt*Leachate_vol_mL/Char_leached_g*1000/1000)

#convert solid elemental concentrations from to g/kg
data_solids<-data_solids %>%
  mutate(Solid_P_g_kg = Solid_P_mg_kg/1000) %>%
  mutate(Solid_Ca_g_kg = Solid_Ca_mg_kg/1000) %>%
  mutate(Solid_Fe_g_kg = Solid_Fe_mg_kg/1000) %>%
  mutate(Solid_Al_g_kg = Solid_Al_mg_kg/1000) %>%
  mutate(Solid_K_g_kg = Solid_K_mg_kg/1000) %>%
  mutate(Solid_Mg_g_kg = Solid_Mg_mg_kg/1000) %>%
  mutate(Solid_Na_g_kg = Solid_Na_mg_kg/1000) %>%
  mutate(Solid_S_g_kg = Solid_S_mg_kg/1000) %>%
  mutate(Solid_C_g_kg = C_weight_per/100*1000) %>%
  mutate(Solid_N_g_kg = N_weight_per/100*1000)

#normalize by P concentration in solid char fix
data<-data %>%
  mutate(Leachate_Ca_perc_solid_part = Leachate_Ca_mg_kgchar_part/Solid_P_mg_kg*100) %>%
  mutate(Leachate_Mg_perc_solid_part = Leachate_Mg_mg_kgchar_part/Solid_P_mg_kg*100) %>%
  mutate(Leachate_P_perc_solid_part = Leachate_P_mg_kgchar_part/Solid_P_mg_kg*100) %>%
  mutate(Leachate_K_perc_solid_part = Leachate_K_mg_kgchar_part/Solid_P_mg_kg*100) %>%
  mutate(Leachate_Na_perc_solid_part = Leachate_Na_mg_kgchar_part/Solid_P_mg_kg*100) %>%
  mutate(Leachate_S_perc_solid_part = Leachate_S_mg_kgchar_part/Solid_P_mg_kg*100) %>%
  mutate(Leachate_Al_perc_solid_part = Leachate_Al_mg_kgchar_part/Solid_P_mg_kg*100) %>%
  mutate(Leachate_Fe_perc_solid_part = Leachate_Fe_mg_kgchar_part/Solid_P_mg_kg*100) %>%
  mutate(Leachate_Ca_perc_solid_filt0.7 = Leachate_Ca_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
  mutate(Leachate_Mg_perc_solid_filt0.7 = Leachate_Mg_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
  mutate(Leachate_P_perc_solid_filt0.7 = Leachate_P_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
  mutate(Leachate_K_perc_solid_filt0.7 = Leachate_K_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
  mutate(Leachate_Na_perc_solid_filt0.7 = Leachate_Na_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
  mutate(Leachate_S_perc_solid_filt0.7 = Leachate_S_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
  mutate(Leachate_Al_perc_solid_filt0.7 = Leachate_Al_mg_kgchar_filt0.7/Solid_P_mg_kg*100) %>%
  mutate(Leachate_Fe_perc_solid_filt0.7 = Leachate_Fe_mg_kgchar_filt0.7/Solid_P_mg_kg*100)

#molar stoichiometry of P:element in solid phase fix
test<-data_solids %>%
  mutate(Fe.P.MolarRatio=(Solid_Fe_mg_kg/55.845/1000)/(Solid_P_mg_kg/30.97/1000))%>%
  mutate(K.P.MolarRatio=(Solid_K_mg_kg/39.0983/1000)/(Solid_P_mg_kg/30.97/1000))%>%
  mutate(Na.P.MolarRatio=(Solid_Na_mg_kg/22.989/1000)/(Solid_P_mg_kg/30.97/1000))%>%
  mutate(Al.P.MolarRatio=(Solid_Al_mg_kg/26.982/1000)/(Solid_P_mg_kg/30.97/1000))%>%
  mutate(Mg.P.MolarRatio=(Solid_Mg_mg_kg/24.305/1000)/(Solid_P_mg_kg/30.97/1000))%>%
  mutate(Ca.P.MolarRatio=(Solid_Ca_mg_kg/40.078/1000)/(Solid_P_mg_kg/30.97/1000))%>%
  mutate(S.P.MolarRatio=(Solid_S_mg_kg/32.065/1000)/(Solid_P_mg_kg/30.97/1000))



#molar stoichiometry of P:element in particulate phase

#molar stoichiometry of P:element in <0.7 phase



#save csv of final collated data with manipulations
write.csv(data, "../data/data.csv", row.names = FALSE) #all data
write.csv(data_solids, "../data/data_solids.csv", row.names = FALSE) #solids only

##Figures and Stats

#import data so you don't have to run above code

data <- read_csv("../data/data.csv")
## Rows: 57 Columns: 99
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (7): Parent_ID, Sample_ID, Land_Coverage_Category, Burn_Treatment, Burn...
## dbl (91): Char_leached_g, Leachate_vol_mL, Burn_Duration, Char_Max_Temp, C_w...
## lgl  (1): FilterSize
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
data_solids <- read_csv("../data/data_solids.csv")
## Rows: 19 Columns: 120
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (10): Parent_ID, Sample_ID, Land_Coverage_Category, Burn_Treatment, Bur...
## dbl (110): Char_leached_g, Leachate_vol_mL, Burn_Duration, Char_Max_Temp, C_...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

####summary stats

#solid elemental concentrations, NMR, and XANES summary
summary_solids_conc <- data_solids %>%
  dplyr::group_by(Burn_Severity, Land_Coverage_Category) %>%
  dplyr::summarize(
    P_Mean = mean(Solid_P_g_kg, na.rm = TRUE),
    P_Standard_Deviation = sd(Solid_P_g_kg, na.rm = TRUE),
    Ca_Mean = mean(Solid_Ca_g_kg, na.rm = TRUE),
    Ca_Standard_Deviation = sd(Solid_Ca_g_kg, na.rm = TRUE),
    Fe_Mean = mean(Solid_Fe_g_kg, na.rm = TRUE),
    Fe_Standard_Deviation = sd(Solid_Fe_g_kg, na.rm = TRUE),
    Al_Mean = mean(Solid_Al_g_kg, na.rm = TRUE),
    Al_Standard_Deviation = sd(Solid_Al_g_kg, na.rm = TRUE),
    K_Mean = mean(Solid_K_g_kg, na.rm = TRUE),
    K_Standard_Deviation = sd(Solid_K_g_kg, na.rm = TRUE),
    Mg_Mean = mean(Solid_Mg_g_kg, na.rm = TRUE),
    Mg_Standard_Deviation = sd(Solid_Mg_g_kg, na.rm = TRUE),
    Na_Mean = mean(Solid_Na_g_kg, na.rm = TRUE),
    Na_Standard_Deviation = sd(Solid_Na_g_kg, na.rm = TRUE),
    S_Mean = mean(Solid_S_g_kg, na.rm = TRUE),
    S_Standard_Deviation = sd(Solid_S_g_kg, na.rm = TRUE),
    C_Mean = mean(Solid_C_g_kg, na.rm = TRUE),
    C_Standard_Deviation = sd(Solid_C_g_kg, na.rm = TRUE),
    N_Mean = mean(Solid_N_g_kg, na.rm = TRUE),
    N_Standard_Deviation = sd(Solid_N_g_kg, na.rm = TRUE),
    Ortho_Mean = mean(Ortho, na.rm = TRUE),
    Ortho_Standard_Deviation = sd(Ortho, na.rm = TRUE),
    Mono_Mean = mean(Mono, na.rm = TRUE),
    Mono_Standard_Deviation = sd(Mono, na.rm = TRUE),
    Di_Mean = mean(Di, na.rm = TRUE),
    Di_Standard_Deviation = sd(Di, na.rm = TRUE),
    Pyro_Mean = mean(Pyro, na.rm = TRUE),
    Pyro_Standard_Deviation = sd(Pyro, na.rm = TRUE),
    XANES_Po_Mean = mean(XANES_Po, na.rm = TRUE),
    XANES_Po_Standard_Deviation = sd(XANES_Po, na.rm = TRUE),
    XANES_Pi_Mean = mean(XANES_Pi, na.rm = TRUE),
    XANES_Pi_Standard_Deviation = sd(XANES_Pi, na.rm = TRUE),
    XANES_Fe_Mean = mean(XANES_Fe, na.rm = TRUE),
    XANES_Fe_Standard_Deviation = sd(XANES_Fe, na.rm = TRUE),
    XANES_Na_Mean = mean(XANES_Na, na.rm = TRUE),
    XANES_Na_Standard_Deviation = sd(XANES_Na, na.rm = TRUE),
    XANES_Ca_Mean = mean(XANES_Ca, na.rm = TRUE),
    XANES_Ca_Standard_Deviation = sd(XANES_Ca, na.rm = TRUE),
    XANES_Al_Mean = mean(XANES_Al, na.rm = TRUE),
    XANES_Al_Standard_Deviation = sd(XANES_Al, na.rm = TRUE),
    XANES_Mg_Mean = mean(XANES_Mg, na.rm = TRUE),
    XANES_Mg_Standard_Deviation = sd(XANES_Mg, na.rm = TRUE),
    XANES_K_Mean = mean(XANES_K, na.rm = TRUE),
    XANES_K_Standard_Deviation = sd(XANES_K, na.rm = TRUE),
    XANES_Pi_Fe_Mean = mean(XANES_Pi_Fe, na.rm = TRUE),
    XANES_Pi_Fe_Standard_Deviation = sd(XANES_Pi_Fe, na.rm = TRUE),
    XANES_Pi_Na_Mean = mean(XANES_Pi_Na, na.rm = TRUE),
    XANES_Pi_Na_Standard_Deviation = sd(XANES_Pi_Na, na.rm = TRUE),
    XANES_Pi_Ca_Mean = mean(XANES_Pi_Ca, na.rm = TRUE),
    XANES_Pi_Ca_Standard_Deviation = sd(XANES_Pi_Ca, na.rm = TRUE),
    XANES_Pi_Al_Mean = mean(XANES_Pi_Al, na.rm = TRUE),
    XANES_Pi_Al_Standard_Deviation = sd(XANES_Pi_Al, na.rm = TRUE),
    XANES_Pi_Mg_Mean = mean(XANES_Pi_Mg, na.rm = TRUE),
    XANES_Pi_Mg_Standard_Deviation = sd(XANES_Pi_Mg, na.rm = TRUE),
    XANES_Pi_K_Mean = mean(XANES_Pi_K, na.rm = TRUE),
    XANES_Pi_K_Standard_Deviation = sd(XANES_Pi_K, na.rm = TRUE),
    XANES_Po_Fe_Mean = mean(XANES_Po_Fe, na.rm = TRUE),
    XANES_Po_Fe_Standard_Deviation = sd(XANES_Po_Fe, na.rm = TRUE),
    XANES_Po_Na_Mean = mean(XANES_Po_Na, na.rm = TRUE),
    XANES_Po_Na_Standard_Deviation = sd(XANES_Po_Na, na.rm = TRUE),
    XANES_Po_Ca_Mean = mean(XANES_Po_Ca, na.rm = TRUE),
    XANES_Po_Ca_Standard_Deviation = sd(XANES_Po_Ca, na.rm = TRUE),
    XANES_Po_Al_Mean = mean(XANES_Po_Al, na.rm = TRUE),
    XANES_Po_Al_Standard_Deviation = sd(XANES_Po_Al, na.rm = TRUE),
    XANES_Po_Mg_Mean = mean(XANES_Po_Mg, na.rm = TRUE),
    XANES_Po_Mg_Standard_Deviation = sd(XANES_Po_Mg, na.rm = TRUE),
    XANES_Po_K_Mean = mean(XANES_Po_K, na.rm = TRUE),
    XANES_Po_K_Standard_Deviation = sd(XANES_Po_K, na.rm = TRUE),
    Count = n()
  )
## `summarise()` has grouped output by 'Burn_Severity'. You can override using the
## `.groups` argument.
write.csv(summary_solids_conc, "../data/summary_solids_conc.csv") #save this for a table


#NMR extraction efficiency
summary_solids <- data_solids %>%
  dplyr::group_by(Burn_Severity, Land_Coverage_Category) %>%
  dplyr::summarize(
    Mean = mean(NMR_EE_P_perc, na.rm = TRUE),
    Standard_Deviation = sd(NMR_EE_P_perc, na.rm = TRUE),
    Count = n()
  )
## `summarise()` has grouped output by 'Burn_Severity'. You can override using the
## `.groups` argument.
#leachate concentrations
P.Conc.leach<- data %>%
  dplyr::group_by(Burn_Severity, Land_Coverage_Category) %>%
  dplyr::summarize(
    Leachate_P_mg_kgchar_part_Mean = mean(Leachate_P_mg_kgchar_part, na.rm = TRUE),
    Leachate_P_mg_kgchar_part_Standard_Deviation = sd(Leachate_P_mg_kgchar_part, na.rm = TRUE),
    Count = n(),
    Leachate_P_mg_kgchar_filt0.7_Mean = mean(Leachate_P_mg_kgchar_filt0.7, na.rm = TRUE),
    Leachate_P_mg_kgchar_filt0.7_Standard_Deviation = sd(Leachate_P_mg_kgchar_filt0.7, na.rm = TRUE),
    Leachate_P_mg_kgchar_unfilt_Mean = mean(Leachate_P_mg_kgchar_unfilt, na.rm = TRUE),
    Leachate_P_mg_kgchar_unfilt_Standard_Deviation = sd(Leachate_P_mg_kgchar_unfilt, na.rm = TRUE)
  )
## `summarise()` has grouped output by 'Burn_Severity'. You can override using the
## `.groups` argument.

####Molybdate

#MBD mg/L figure
ggplot(data, aes(x=Burn_Severity, y=MBD_P_mg_per_L_filt0.7, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Soluble~Reactive~P~(mg~P~L^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))
## Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).

#MBD mg/kg figure; likely use this one for publication
MBD.BS.fig<-ggplot(data, aes(x=Burn_Severity, y=MBD_P_mg_kgchar/1000, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Soluble~Reactive~P~(g~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))

cowplot::save_plot("../figures/MBD.BS.fig.pdf", MBD.BS.fig, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)
## Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).
#MBD as percent of solid P concentration figure
ggplot(data, aes(x=Burn_Severity, y=MBD_P_percentofsolid, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab("Soluble Reactive P (%)")+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))
## Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).

#MBD g/kg stats; likely use this one for publication
model<-lmer(log10(MBD_P_mg_kgchar/1000)~Burn_Severity*Land_Coverage_Category+(1|Parent_ID), data=data) #if want to use maximum likelihood , REML=FALSE
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
model.BS<-lmer(log10(MBD_P_mg_kgchar/1000)~Land_Coverage_Category+(1|Parent_ID), data=data)
model.LC<-lmer(log10(MBD_P_mg_kgchar/1000)~Burn_Severity+(1|Parent_ID), data=data)
model.int<-lmer(log10(MBD_P_mg_kgchar/1000)~Burn_Severity+Land_Coverage_Category+(1|Parent_ID), data=data)
anova(model,model.BS) #not sig
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.BS: log10(MBD_P_mg_kgchar/1000) ~ Land_Coverage_Category + (1 | Parent_ID)
## model: log10(MBD_P_mg_kgchar/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##          npar    AIC    BIC   logLik deviance Chisq Df Pr(>Chisq)
## model.BS    4 30.181 38.137 -11.0905   22.181                    
## model       9 32.325 50.226  -7.1625   14.325 7.856  5     0.1644
anova(model,model.LC) #sig
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.LC: log10(MBD_P_mg_kgchar/1000) ~ Burn_Severity + (1 | Parent_ID)
## model: log10(MBD_P_mg_kgchar/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##          npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)  
## model.LC    6 37.291 49.225 -12.6453   25.291                       
## model       9 32.325 50.226  -7.1625   14.325 10.966  3    0.01191 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model,model.int) #sig
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.int: log10(MBD_P_mg_kgchar/1000) ~ Burn_Severity + Land_Coverage_Category + (1 | Parent_ID)
## model: log10(MBD_P_mg_kgchar/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##           npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)  
## model.int    7 35.410 49.333 -10.7051   21.410                       
## model        9 32.325 50.226  -7.1625   14.325 7.0852  2    0.02894 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#check VIF
vif(model)
##                                          GVIF Df GVIF^(1/(2*Df))
## Burn_Severity                        2.567759  3        1.170197
## Land_Coverage_Category               5.502341  1        2.345707
## Burn_Severity:Land_Coverage_Category 9.075768  2        1.735685
#model assumptions
hist(log10(data$MBD_P_mg_kgchar))

#linearity
ggqqplot(residuals(model))

#check normality of residuals 
shapiro.test(log10(data$MBD_P_mg_kgchar)) #p > 0.05 so can assume normality 
## 
##  Shapiro-Wilk normality test
## 
## data:  log10(data$MBD_P_mg_kgchar)
## W = 0.95892, p-value = 0.06167
#check homogeneity of variance assumption (homoscedastcity)
plot(fitted(model),residuals(model))

performance::check_heteroscedasticity(model) #p > 0.05 so can assume homoscedastic
## OK: Error variance appears to be homoscedastic (p = 0.632).
#post hoc tests for interaction term
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Unburned"),adjust="tukey")
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE   df lower.CL upper.CL
##  Unburned      Douglas-fir forest     -1.183 0.221 11.7    -1.67  -0.6994
##  Unburned      Sagebrush shrubland    -0.701 0.313 11.7    -1.39  -0.0163
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                     estimate    SE
##  (Unburned Douglas-fir forest) - Unburned Sagebrush shrubland   -0.483 0.384
##    df t.ratio p.value
##  11.7  -1.258  0.2328
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale.
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Low"),adjust="tukey")
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE   df lower.CL upper.CL
##  Low           Douglas-fir forest     -1.266 0.140 11.7    -1.57   -0.960
##  Low           Sagebrush shrubland    -0.551 0.221 11.7    -1.03   -0.067
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                           estimate    SE   df t.ratio
##  (Low Douglas-fir forest) - Low Sagebrush shrubland   -0.715 0.262 11.7  -2.729
##  p.value
##   0.0187
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale.
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Moderate"),adjust="tukey")
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE   df lower.CL upper.CL
##  Moderate      Douglas-fir forest      -1.05 0.181 11.7    -1.45   -0.655
##  Moderate      Sagebrush shrubland     -1.23 0.221 11.7    -1.71   -0.742
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                     estimate    SE
##  (Moderate Douglas-fir forest) - Moderate Sagebrush shrubland    0.176 0.286
##    df t.ratio p.value
##  11.7   0.616  0.5500
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale.
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Douglas-fir forest"),adjust="tukey")
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE   df lower.CL upper.CL
##  High          Douglas-fir forest      -1.07 0.161 13.1    -1.42   -0.726
##  Low           Douglas-fir forest      -1.27 0.140 11.7    -1.57   -0.960
##  Moderate      Douglas-fir forest      -1.05 0.181 11.7    -1.45   -0.655
##  Unburned      Douglas-fir forest      -1.18 0.221 11.7    -1.67   -0.699
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                      estimate    SE
##  (High Douglas-fir forest) - (Low Douglas-fir forest)            0.1918 0.214
##  (High Douglas-fir forest) - (Moderate Douglas-fir forest)      -0.0238 0.242
##  (High Douglas-fir forest) - (Unburned Douglas-fir forest)       0.1093 0.274
##  (Low Douglas-fir forest) - (Moderate Douglas-fir forest)       -0.2156 0.229
##  (Low Douglas-fir forest) - (Unburned Douglas-fir forest)       -0.0825 0.262
##  (Moderate Douglas-fir forest) - (Unburned Douglas-fir forest)   0.1331 0.286
##    df t.ratio p.value
##  12.4   0.898  0.8062
##  12.3  -0.098  0.9996
##  12.1   0.399  0.9775
##  11.7  -0.943  0.7830
##  11.7  -0.315  0.9886
##  11.7   0.466  0.9651
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Sagebrush shrubland"),adjust="tukey")
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE   df lower.CL upper.CL
##  High          Sagebrush shrubland    nonEst    NA   NA       NA       NA
##  Low           Sagebrush shrubland    -0.551 0.221 11.7    -1.03  -0.0670
##  Moderate      Sagebrush shrubland    -1.226 0.221 11.7    -1.71  -0.7423
##  Unburned      Sagebrush shrubland    -0.701 0.313 11.7    -1.39  -0.0163
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                    estimate    SE
##  High Sagebrush shrubland - Low Sagebrush shrubland            nonEst    NA
##  High Sagebrush shrubland - Moderate Sagebrush shrubland       nonEst    NA
##  High Sagebrush shrubland - Unburned Sagebrush shrubland       nonEst    NA
##  Low Sagebrush shrubland - Moderate Sagebrush shrubland         0.675 0.313
##  Low Sagebrush shrubland - Unburned Sagebrush shrubland         0.150 0.384
##  Moderate Sagebrush shrubland - Unburned Sagebrush shrubland   -0.526 0.384
##    df t.ratio p.value
##    NA      NA      NA
##    NA      NA      NA
##    NA      NA      NA
##  11.7   2.157  0.1908
##  11.7   0.391  0.9788
##  11.7  -1.370  0.5396
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates

####ICP Leachates

#convert to percent of total P in unfiltered sample
#ICP$P.filt.percentofunfilt<-(ICP$Total_P_mg_per_L_filt/ICP$Total_P_mg_per_L_unfilt)*100
#ICP$P.part.percentofunfilt<-(ICP$Total_P_mg_per_L_part/ICP$Total_P_mg_per_L_unfilt)*100

#total P in unfiltered leachate
Leachate.total.P.BS.boxplot <- ggplot(data, aes(x=Burn_Severity, y=Leachate_P_mg_kgchar_unfilt, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Total~P~(mg~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))

cowplot::save_plot("../figures/Leachate.total.P.BS.boxplot.pdf", Leachate.total.P.BS.boxplot, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)

Leachate.filt0.7.total.P.BS.boxplot <- ggplot(data, aes(x=Burn_Severity, y=Leachate_P_mg_kgchar_filt0.7, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Total~Soluble~P~(mg~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))

cowplot::save_plot("../figures/Leachate.filt0.7.total.P.BS.boxplot.pdf", Leachate.filt0.7.total.P.BS.boxplot, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)

Leachate.unfilt.total.P.BS.boxplot <- ggplot(data, aes(x=Burn_Severity, y=Leachate_P_mg_kgchar_part, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Total~Particulate~P~(mg~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))

cowplot::save_plot("../figures/Leachate.unfilt.total.P.BS.boxplot.pdf", Leachate.unfilt.total.P.BS.boxplot, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)


#stacked bar chart
P.Conc.leach.bar <- P.Conc.leach %>%
  mutate(sd.high=Leachate_P_mg_kgchar_unfilt_Mean+Leachate_P_mg_kgchar_unfilt_Standard_Deviation) %>%
  mutate(sd.low=Leachate_P_mg_kgchar_unfilt_Mean-Leachate_P_mg_kgchar_unfilt_Standard_Deviation)


P.Conc.leach_melt<-melt(P.Conc.leach.bar, id=c("Burn_Severity", "Land_Coverage_Category","sd.low", "sd.high")) %>%
  filter(variable %in% c("Leachate_P_mg_kgchar_part_Mean" , "Leachate_P_mg_kgchar_filt0.7_Mean"))


Leachate_stacked_bar <- ggplot(subset(P.Conc.leach_melt, variable=="Leachate_P_mg_kgchar_part_Mean" | variable=="Leachate_P_mg_kgchar_filt0.7_Mean"), aes(x = Burn_Severity, y = value, fill=variable)) + geom_bar(stat = "identity")  + theme_classic()+ylab(expression(paste(Leachate~Mean~Total~P~(mg~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#2274A5", "#E7DFC6"), labels = c("Particulate", "Soluble"), name = "")+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15)) + geom_errorbar(aes(ymax=sd.high, ymin=sd.low), width = 0.2, size = 1, color = "black")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
cowplot::save_plot("../figures/Leachate_stacked_bar.pdf", Leachate_stacked_bar, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)


#ICP unfiltered total P g/kg stats; likely use this one for publication
model<-lmer(log10(Leachate_P_mg_kgchar_unfilt/1000)~Burn_Severity*Land_Coverage_Category+(1|Parent_ID), data=data)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
model.BS<-lmer(log10(Leachate_P_mg_kgchar_unfilt/1000)~Land_Coverage_Category+(1|Parent_ID), data=data)
model.LC<-lmer(log10(Leachate_P_mg_kgchar_unfilt/1000)~Burn_Severity+(1|Parent_ID), data=data)
model.int<-lmer(log10(Leachate_P_mg_kgchar_unfilt/1000)~Burn_Severity+Land_Coverage_Category+(1|Parent_ID), data=data)
anova(model,model.BS)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.BS: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##          npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## model.BS    4 44.894 53.066 -18.4470   36.894                         
## model       9 24.615 43.003  -3.3077    6.615 30.279  5    1.3e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model,model.LC)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.LC: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Burn_Severity + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##          npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## model.LC    6 40.336 52.595 -14.1681  28.3363                         
## model       9 24.615 43.003  -3.3077   6.6153 21.721  3  7.456e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model,model.int)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.int: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Burn_Severity + Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_unfilt/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##           npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## model.int    7 26.705 41.006 -6.3524  12.7048                       
## model        9 24.615 43.003 -3.3077   6.6153 6.0895  2    0.04761 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#check VIF
vif(model)
##                                          GVIF Df GVIF^(1/(2*Df))
## Burn_Severity                        2.578947  3        1.171045
## Land_Coverage_Category               5.526316  1        2.350812
## Burn_Severity:Land_Coverage_Category 9.105263  2        1.737093
#model assumptions
hist(log10(data$Leachate_P_mg_kgchar_unfilt/1000))

#linearity
ggqqplot(residuals(model))

#check normality of residuals 
shapiro.test(log10(data$Leachate_P_mg_kgchar_unfilt)) #p > 0.05 so can assume normality 
## 
##  Shapiro-Wilk normality test
## 
## data:  log10(data$Leachate_P_mg_kgchar_unfilt)
## W = 0.96458, p-value = 0.09345
#check homogeneity of variance assumption (homoscedastcity)
plot(fitted(model),residuals(model))

performance::check_heteroscedasticity(model) #p > 0.05 so can assume homoscedastic
## OK: Error variance appears to be homoscedastic (p = 0.913).
#post hoc tests for interaction term
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Unburned"),adjust="tukey") 
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE df lower.CL upper.CL
##  Unburned      Douglas-fir forest     -0.885 0.231 12    -1.39   -0.381
##  Unburned      Sagebrush shrubland    -0.501 0.327 12    -1.21    0.211
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                     estimate  SE df
##  (Unburned Douglas-fir forest) - Unburned Sagebrush shrubland   -0.384 0.4 12
##  t.ratio p.value
##   -0.958  0.3571
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale.
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Low"),adjust="tukey")
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE df lower.CL upper.CL
##  Low           Douglas-fir forest     -0.843 0.146 12   -1.162   -0.524
##  Low           Sagebrush shrubland    -0.223 0.231 12   -0.726    0.281
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                           estimate    SE df t.ratio
##  (Low Douglas-fir forest) - Low Sagebrush shrubland    -0.62 0.274 12  -2.268
##  p.value
##   0.0426
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale.
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Moderate"),adjust="tukey")
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE df lower.CL upper.CL
##  Moderate      Douglas-fir forest     -0.488 0.189 12   -0.899  -0.0765
##  Moderate      Sagebrush shrubland     0.833 0.231 12    0.329   1.3368
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                     estimate    SE df
##  (Moderate Douglas-fir forest) - Moderate Sagebrush shrubland    -1.32 0.298 12
##  t.ratio p.value
##   -4.425  0.0008
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale.
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Douglas-fir forest"),adjust="tukey")
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE df lower.CL upper.CL
##  High          Douglas-fir forest      0.299 0.163 12  -0.0575   0.6549
##  Low           Douglas-fir forest     -0.843 0.146 12  -1.1616  -0.5244
##  Moderate      Douglas-fir forest     -0.488 0.189 12  -0.8991  -0.0765
##  Unburned      Douglas-fir forest     -0.885 0.231 12  -1.3885  -0.3810
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                      estimate    SE
##  (High Douglas-fir forest) - (Low Douglas-fir forest)            1.1417 0.219
##  (High Douglas-fir forest) - (Moderate Douglas-fir forest)       0.7864 0.250
##  (High Douglas-fir forest) - (Unburned Douglas-fir forest)       1.1834 0.283
##  (Low Douglas-fir forest) - (Moderate Douglas-fir forest)       -0.3552 0.239
##  (Low Douglas-fir forest) - (Unburned Douglas-fir forest)        0.0418 0.274
##  (Moderate Douglas-fir forest) - (Unburned Douglas-fir forest)   0.3970 0.298
##  df t.ratio p.value
##  12   5.205  0.0011
##  12   3.149  0.0366
##  12   4.179  0.0061
##  12  -1.488  0.4739
##  12   0.153  0.9987
##  12   1.330  0.5627
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Sagebrush shrubland"),adjust="tukey")
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE df lower.CL upper.CL
##  High          Sagebrush shrubland    nonEst    NA NA       NA       NA
##  Low           Sagebrush shrubland    -0.223 0.231 12   -0.726    0.281
##  Moderate      Sagebrush shrubland     0.833 0.231 12    0.329    1.337
##  Unburned      Sagebrush shrubland    -0.501 0.327 12   -1.214    0.211
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                    estimate    SE df
##  High Sagebrush shrubland - Low Sagebrush shrubland            nonEst    NA NA
##  High Sagebrush shrubland - Moderate Sagebrush shrubland       nonEst    NA NA
##  High Sagebrush shrubland - Unburned Sagebrush shrubland       nonEst    NA NA
##  Low Sagebrush shrubland - Moderate Sagebrush shrubland        -1.056 0.327 12
##  Low Sagebrush shrubland - Unburned Sagebrush shrubland         0.279 0.400 12
##  Moderate Sagebrush shrubland - Unburned Sagebrush shrubland    1.334 0.400 12
##  t.ratio p.value
##       NA      NA
##       NA      NA
##       NA      NA
##   -3.228  0.0318
##    0.696  0.8967
##    3.332  0.0266
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates
#ICP filtered (<0.7 nominal phase) total P stats; likely use this one for publication
model<-lmer(log10(Leachate_P_mg_kgchar_filt0.7/1000)~Burn_Severity*Land_Coverage_Category+(1|Parent_ID), data=data)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
model.BS<-lmer(log10(Leachate_P_mg_kgchar_filt0.7/1000)~Land_Coverage_Category+(1|Parent_ID), data=data)
model.LC<-lmer(log10(Leachate_P_mg_kgchar_filt0.7/1000)~Burn_Severity+(1|Parent_ID), data=data)
model.int<-lmer(log10(Leachate_P_mg_kgchar_filt0.7/1000)~Burn_Severity+Land_Coverage_Category+(1|Parent_ID), data=data)
anova(model,model.BS) #not sig
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.BS: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##          npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## model.BS    4 -37.826 -29.654 22.913  -45.826                     
## model       9 -33.886 -15.498 25.943  -51.886 6.0598  5     0.3004
anova(model,model.LC) #sig
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.LC: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Burn_Severity + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##          npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)   
## model.LC    6 -26.659 -14.401 19.330  -38.659                        
## model       9 -33.886 -15.498 25.943  -51.886 13.226  3   0.004172 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model,model.int) #not sig but since BS also isn't sig I probably don't care that much about the LC main effect (so not being able to run the post hoc isn't an issue)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.int: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Burn_Severity + Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_filt0.7/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##           npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)  
## model.int    7 -32.512 -18.211 23.256  -46.512                       
## model        9 -33.886 -15.498 25.943  -51.886 5.3737  2    0.06809 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#check VIF
vif(model)
##                                          GVIF Df GVIF^(1/(2*Df))
## Burn_Severity                        2.578947  3        1.171045
## Land_Coverage_Category               5.526316  1        2.350812
## Burn_Severity:Land_Coverage_Category 9.105263  2        1.737093
#model assumptions
hist(log10(data$Leachate_P_mg_kgchar_filt0.7/1000))

ggqqplot(residuals(model))

#check normality of residuals 
shapiro.test(log10(data$Leachate_P_mg_kgchar_unfilt)) #p > 0.05 so can assume normality 
## 
##  Shapiro-Wilk normality test
## 
## data:  log10(data$Leachate_P_mg_kgchar_unfilt)
## W = 0.96458, p-value = 0.09345
#check homogeneity of variance assumption (homoscedastcity)
plot(fitted(model),residuals(model))

performance::check_heteroscedasticity(model) #need to look into this (fix)
## Warning in sqrt(insight::get_deviance(x)/(insight::n_obs(x) -
## sum(!is.na(estimates)))): NaNs produced
## Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).
#ICP particulate (>0.7 nominal phase) total P stats; likely use this one for publication
model<-lmer(log10(Leachate_P_mg_kgchar_part/1000)~Burn_Severity*Land_Coverage_Category+(1|Parent_ID), data=data)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
model.BS<-lmer(log10(Leachate_P_mg_kgchar_part/1000)~Land_Coverage_Category+(1|Parent_ID), data=data)
model.LC<-lmer(log10(Leachate_P_mg_kgchar_part/1000)~Burn_Severity+(1|Parent_ID), data=data)
model.int<-lmer(log10(Leachate_P_mg_kgchar_part/1000)~Burn_Severity+Land_Coverage_Category+(1|Parent_ID), data=data)
anova(model,model.BS) #sig
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.BS: log10(Leachate_P_mg_kgchar_part/1000) ~ Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_part/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## model.BS    4 116.60 124.77 -54.298   108.60                         
## model       9  90.37 108.76 -36.185    72.37 36.226  5  8.559e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model,model.LC) #sig
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.LC: log10(Leachate_P_mg_kgchar_part/1000) ~ Burn_Severity + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_part/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## model.LC    6 101.81 114.07 -44.905   89.811                         
## model       9  90.37 108.76 -36.185   72.370 17.441  3  0.0005736 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model,model.int) #sig
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## model.int: log10(Leachate_P_mg_kgchar_part/1000) ~ Burn_Severity + Land_Coverage_Category + (1 | Parent_ID)
## model: log10(Leachate_P_mg_kgchar_part/1000) ~ Burn_Severity * Land_Coverage_Category + (1 | Parent_ID)
##           npar   AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)   
## model.int    7 98.75 113.05 -42.375    84.75                       
## model        9 90.37 108.76 -36.185    72.37 12.38  2    0.00205 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#check VIF
vif(model)
##                                          GVIF Df GVIF^(1/(2*Df))
## Burn_Severity                        2.578947  3        1.171045
## Land_Coverage_Category               5.526316  1        2.350812
## Burn_Severity:Land_Coverage_Category 9.105263  2        1.737093
#model assumptions
hist(log10(data$Leachate_P_mg_kgchar_filt0.7/1000))

ggqqplot(residuals(model))

#check normality of residuals 
shapiro.test(log10(data$Leachate_P_mg_kgchar_unfilt)) #p > 0.05 so can assume normality 
## 
##  Shapiro-Wilk normality test
## 
## data:  log10(data$Leachate_P_mg_kgchar_unfilt)
## W = 0.96458, p-value = 0.09345
#check homogeneity of variance assumption (homoscedastcity)
plot(fitted(model),residuals(model))

performance::check_heteroscedasticity(model) #p > 0.05 so can assume homoscedastic
## OK: Error variance appears to be homoscedastic (p = 0.745).
#interaction post hoc: 
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Unburned"),adjust="tukey") 
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE df lower.CL upper.CL
##  Unburned      Douglas-fir forest      -1.46 0.310 12    -2.14   -0.788
##  Unburned      Sagebrush shrubland     -2.16 0.438 12    -3.11   -1.203
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                     estimate    SE df
##  (Unburned Douglas-fir forest) - Unburned Sagebrush shrubland    0.696 0.537 12
##  t.ratio p.value
##    1.296  0.2195
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale.
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Low"),adjust="tukey") #between 
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE df lower.CL upper.CL
##  Low           Douglas-fir forest     -1.322 0.196 12    -1.75   -0.895
##  Low           Sagebrush shrubland    -0.791 0.310 12    -1.47   -0.115
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                           estimate    SE df t.ratio
##  (Low Douglas-fir forest) - Low Sagebrush shrubland   -0.531 0.367 12  -1.449
##  p.value
##   0.1731
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale.
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Moderate"),adjust="tukey") 
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE df lower.CL upper.CL
##  Moderate      Douglas-fir forest     -0.687 0.253 12   -1.238   -0.136
##  Moderate      Sagebrush shrubland     0.825 0.310 12    0.149    1.500
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                     estimate  SE df
##  (Moderate Douglas-fir forest) - Moderate Sagebrush shrubland    -1.51 0.4 12
##  t.ratio p.value
##   -3.778  0.0026
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale.
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Douglas-fir forest"),adjust="tukey") 
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE df lower.CL upper.CL
##  High          Douglas-fir forest      0.269 0.219 12   -0.208    0.747
##  Low           Douglas-fir forest     -1.322 0.196 12   -1.749   -0.895
##  Moderate      Douglas-fir forest     -0.687 0.253 12   -1.238   -0.136
##  Unburned      Douglas-fir forest     -1.463 0.310 12   -2.138   -0.788
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                      estimate    SE
##  (High Douglas-fir forest) - (Low Douglas-fir forest)             1.591 0.294
##  (High Douglas-fir forest) - (Moderate Douglas-fir forest)        0.956 0.335
##  (High Douglas-fir forest) - (Unburned Douglas-fir forest)        1.732 0.380
##  (Low Douglas-fir forest) - (Moderate Douglas-fir forest)        -0.635 0.320
##  (Low Douglas-fir forest) - (Unburned Douglas-fir forest)         0.141 0.367
##  (Moderate Douglas-fir forest) - (Unburned Douglas-fir forest)    0.776 0.400
##  df t.ratio p.value
##  12   5.410  0.0008
##  12   2.856  0.0608
##  12   4.563  0.0031
##  12  -1.983  0.2472
##  12   0.385  0.9797
##  12   1.939  0.2633
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates
emmeans::lsmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Sagebrush shrubland"),adjust="tukey") 
## $lsmeans
##  Burn_Severity Land_Coverage_Category lsmean    SE df lower.CL upper.CL
##  High          Sagebrush shrubland    nonEst    NA NA       NA       NA
##  Low           Sagebrush shrubland    -0.791 0.310 12   -1.466   -0.115
##  Moderate      Sagebrush shrubland     0.825 0.310 12    0.149    1.500
##  Unburned      Sagebrush shrubland    -2.159 0.438 12   -3.114   -1.203
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                    estimate    SE df
##  High Sagebrush shrubland - Low Sagebrush shrubland            nonEst    NA NA
##  High Sagebrush shrubland - Moderate Sagebrush shrubland       nonEst    NA NA
##  High Sagebrush shrubland - Unburned Sagebrush shrubland       nonEst    NA NA
##  Low Sagebrush shrubland - Moderate Sagebrush shrubland         -1.62 0.438 12
##  Low Sagebrush shrubland - Unburned Sagebrush shrubland          1.37 0.537 12
##  Moderate Sagebrush shrubland - Unburned Sagebrush shrubland     2.98 0.537 12
##  t.ratio p.value
##       NA      NA
##       NA      NA
##       NA      NA
##   -3.685  0.0143
##    2.548  0.1020
##    5.557  0.0006
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates

####Solid ICP

solid.P.BS.fig<-ggplot(data_solids, aes(x=Burn_Severity, y=Solid_P_mg_kg, fill=Land_Coverage_Category))+geom_boxplot(show.legend=FALSE)+theme_classic()+ylab(expression(paste(Solid~Total~P~(mg~P~kg^-1), sep="")))+ xlab("Burn Severity")+scale_fill_manual(values=c("#346648", "#dbb40d"))+facet_grid(.~Land_Coverage_Category)+theme(text = element_text(size = 15))

cowplot::save_plot("../figures/solid.P.BS.fig.pdf", solid.P.BS.fig, ncol = 1, nrow = 1, base_aspect_ratio= 2:1, dpi=72)

#ICP solid total P stats; likely use this one for publication
model<-aov(log10(Solid_P_mg_kg)~Burn_Severity*Land_Coverage_Category, data = data_solids)
summary(model)
##                                      Df Sum Sq Mean Sq F value   Pr(>F)    
## Burn_Severity                         3 1.0646  0.3549  16.124 0.000165 ***
## Land_Coverage_Category                1 0.7795  0.7795  35.419  6.7e-05 ***
## Burn_Severity:Land_Coverage_Category  2 0.2741  0.1370   6.227 0.013965 *  
## Residuals                            12 0.2641  0.0220                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#model assumptions
hist(log10(data_solids$Solid_P_mg_kg))

ggqqplot(residuals(model))

#check normality of residuals 
shapiro.test(log10(data_solids$Solid_P_mg_kg)) #p > 0.05 so can assume normality 
## 
##  Shapiro-Wilk normality test
## 
## data:  log10(data_solids$Solid_P_mg_kg)
## W = 0.94907, p-value = 0.3811
#check homogeneity of variance assumption (homoscedastcity)
plot(fitted(model),residuals(model))

performance::check_heteroscedasticity(model) #p > 0.05 so can assume homoscedastic
## OK: Error variance appears to be homoscedastic (p = 0.558).
#interaction post hoc: 
emmeans::emmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Unburned"),adjust="tukey") 
## $emmeans
##  Burn_Severity Land_Coverage_Category emmean    SE df lower.CL upper.CL
##  Unburned      Douglas-fir forest       3.09 0.105 12     2.86     3.32
##  Unburned      Sagebrush shrubland      3.10 0.148 12     2.78     3.43
## 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                     estimate    SE df
##  (Unburned Douglas-fir forest) - Unburned Sagebrush shrubland -0.00881 0.182 12
##  t.ratio p.value
##   -0.048  0.9621
## 
## Results are given on the log10 (not the response) scale.
emmeans::emmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Low"),adjust="tukey") #between 
## $emmeans
##  Burn_Severity Land_Coverage_Category emmean     SE df lower.CL upper.CL
##  Low           Douglas-fir forest       3.24 0.0663 12     3.09     3.38
##  Low           Sagebrush shrubland      3.68 0.1049 12     3.45     3.91
## 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                           estimate    SE df t.ratio
##  (Low Douglas-fir forest) - Low Sagebrush shrubland   -0.444 0.124 12  -3.579
##  p.value
##   0.0038
## 
## Results are given on the log10 (not the response) scale.
emmeans::emmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Burn_Severity="Moderate"),adjust="tukey") 
## $emmeans
##  Burn_Severity Land_Coverage_Category emmean     SE df lower.CL upper.CL
##  Moderate      Douglas-fir forest       3.35 0.0857 12     3.16     3.54
##  Moderate      Sagebrush shrubland      4.15 0.1049 12     3.92     4.38
## 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                     estimate    SE df
##  (Moderate Douglas-fir forest) - Moderate Sagebrush shrubland   -0.802 0.135 12
##  t.ratio p.value
##   -5.922  0.0001
## 
## Results are given on the log10 (not the response) scale.
emmeans::emmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Douglas-fir forest"),adjust="tukey") 
## $emmeans
##  Burn_Severity Land_Coverage_Category emmean     SE df lower.CL upper.CL
##  High          Douglas-fir forest       3.78 0.0742 12     3.61     3.94
##  Low           Douglas-fir forest       3.24 0.0663 12     3.09     3.38
##  Moderate      Douglas-fir forest       3.35 0.0857 12     3.16     3.54
##  Unburned      Douglas-fir forest       3.09 0.1049 12     2.86     3.32
## 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                      estimate     SE
##  (High Douglas-fir forest) - (Low Douglas-fir forest)             0.537 0.0995
##  (High Douglas-fir forest) - (Moderate Douglas-fir forest)        0.424 0.1133
##  (High Douglas-fir forest) - (Unburned Douglas-fir forest)        0.682 0.1285
##  (Low Douglas-fir forest) - (Moderate Douglas-fir forest)        -0.113 0.1083
##  (Low Douglas-fir forest) - (Unburned Douglas-fir forest)         0.145 0.1241
##  (Moderate Douglas-fir forest) - (Unburned Douglas-fir forest)    0.258 0.1354
##  df t.ratio p.value
##  12   5.394  0.0008
##  12   3.739  0.0130
##  12   5.305  0.0009
##  12  -1.044  0.7281
##  12   1.167  0.6578
##  12   1.904  0.2768
## 
## Results are given on the log10 (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates
emmeans::emmeans(model, pairwise~Burn_Severity*Land_Coverage_Category, at=list(Land_Coverage_Category="Sagebrush shrubland"),adjust="tukey") 
## $emmeans
##  Burn_Severity Land_Coverage_Category emmean    SE df lower.CL upper.CL
##  High          Sagebrush shrubland    nonEst    NA NA       NA       NA
##  Low           Sagebrush shrubland      3.68 0.105 12     3.45     3.91
##  Moderate      Sagebrush shrubland      4.15 0.105 12     3.92     4.38
##  Unburned      Sagebrush shrubland      3.10 0.148 12     2.78     3.43
## 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                    estimate    SE df
##  High Sagebrush shrubland - Low Sagebrush shrubland            nonEst    NA NA
##  High Sagebrush shrubland - Moderate Sagebrush shrubland       nonEst    NA NA
##  High Sagebrush shrubland - Unburned Sagebrush shrubland       nonEst    NA NA
##  Low Sagebrush shrubland - Moderate Sagebrush shrubland        -0.471 0.148 12
##  Low Sagebrush shrubland - Unburned Sagebrush shrubland         0.580 0.182 12
##  Moderate Sagebrush shrubland - Unburned Sagebrush shrubland    1.051 0.182 12
##  t.ratio p.value
##       NA      NA
##       NA      NA
##       NA      NA
##   -3.174  0.0350
##    3.193  0.0339
##    5.784  0.0004
## 
## Results are given on the log10 (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates

####Correlations

#across both feedstocks
#subset only solid elemental concentrations (in mg/kg)
data_solids_corr <- data_solids %>%
  select(matches("Solid"))

# Rename columns by extracting everything between "Solid_" and "mg_Kg"
data_solids_corr <- data_solids_corr %>%
  rename_with(
    ~sub("Solid_(.*?)_mg_kg", "\\1", .),
    starts_with("Solid_")) %>%
  select(P, C, N, Ca, Mg, Na, Fe, Al, S, K)

#comput r=the correlation matrix, n=the matrix of the number of observations used in analyzing each pari fo variables, and P=the p-values corresponding to the significance levels of correlations
cor_5 <- rcorr(as.matrix(data_solids_corr), type=c("spearman"))
#define M as the correlation coefficient
M <- cor_5$r
#define p_mat as the p-value
p_mat <- cor_5$P 
#plot only significant pearson correlation coeffcients when p value >0.05
corrplot(M, type = "upper", p.mat = p_mat, sig.level = 0.05, insig = "blank", tl.col = "black", diag=FALSE, na.label="NA", tl.cex=1, cl.cex = 1, number.cex = 0.8, addCoef.col = "white", col = colorRampPalette(c("#346648", "#dbb40d"))(100))

#Douglas-fir forest
#subset only solid elemental concentrations (in mg/kg)
data_solids_corr <- data_solids  %>%
  filter(str_detect(Land_Coverage_Category, "Douglas")) %>%
  select(matches("Solid"))

# Rename columns by extracting everything between "Solid_" and "mg_Kg"
data_solids_corr <- data_solids_corr %>%
  rename_with(
    ~sub("Solid_(.*?)_mg_kg", "\\1", .),
    starts_with("Solid_")) %>%
  select(P, C, N, Ca, Mg, Na, Fe, Al, S, K)

#comput r=the correlation matrix, n=the matrix of the number of observations used in analyzing each pari fo variables, and P=the p-values corresponding to the significance levels of correlations
cor_5 <- rcorr(as.matrix(data_solids_corr), type=c("spearman"))
#define M as the correlation coefficient
M <- cor_5$r
#define p_mat as the p-value
p_mat <- cor_5$P 
#plot only significant pearson correlation coeffcients when p value >0.05
corrplot(M, type = "upper", p.mat = p_mat, sig.level = 0.05, insig = "blank", tl.col = "black", diag=FALSE, na.label="NA", tl.cex=1, cl.cex = 1, number.cex = 0.8, addCoef.col = "white", col = colorRampPalette(c("#346648", "#dbb40d"))(100))

#Big sagebrush
#subset only solid elemental concentrations (in mg/kg)
data_solids_corr <- data_solids  %>%
  filter(str_detect(Land_Coverage_Category, "Sagebrush")) %>%
  select(matches("Solid"))

# Rename columns by extracting everything between "Solid_" and "mg_Kg"
data_solids_corr <- data_solids_corr %>%
  rename_with(
    ~sub("Solid_(.*?)_mg_kg", "\\1", .),
    starts_with("Solid_")) %>%
  select(P, C, N, Ca, Mg, Na, Fe, Al, S, K)

#comput r=the correlation matrix, n=the matrix of the number of observations used in analyzing each pari fo variables, and P=the p-values corresponding to the significance levels of correlations
cor_5 <- rcorr(as.matrix(data_solids_corr), type=c("spearman"))
#define M as the correlation coefficient
M <- cor_5$r
#define p_mat as the p-value
p_mat <- cor_5$P 
#plot only significant pearson correlation coeffcients when p value >0.05
corrplot(M, type = "upper", p.mat = p_mat, sig.level = 0.05, insig = "blank", tl.col = "black", diag=FALSE, na.label="NA", tl.cex=1, cl.cex = 1, number.cex = 0.8, addCoef.col = "white", col = colorRampPalette(c("#346648", "#dbb40d"))(100))

#P and C
ggplot(data_solids, aes(x=Solid_C_mg_kg, y=Solid_P_mg_kg, color = Land_Coverage_Category)) + geom_point() + geom_smooth(method="lm", se = FALSE) + theme_classic() + scale_color_manual(values = c("#346648", "#dbb40d"))
## `geom_smooth()` using formula = 'y ~ x'

#all feedstocks
model<-lm(log(Solid_C_mg_kg)~log(Solid_P_mg_kg), data = data_solids)
summary(model)
## 
## Call:
## lm(formula = log(Solid_C_mg_kg) ~ log(Solid_P_mg_kg), data = data_solids)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1667 -0.3490  0.1323  0.3302  0.5910 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         13.5156     1.0777   12.54 5.11e-10 ***
## log(Solid_P_mg_kg)  -0.1815     0.1334   -1.36    0.191    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4742 on 17 degrees of freedom
## Multiple R-squared:  0.09818,    Adjusted R-squared:  0.04513 
## F-statistic: 1.851 on 1 and 17 DF,  p-value: 0.1915
#Douglas fir only
model<-lm(log(Solid_C_mg_kg)~log(Solid_P_mg_kg), data = subset(data_solids, Land_Coverage_Category=="Douglas-fir forest"))
summary(model)
## 
## Call:
## lm(formula = log(Solid_C_mg_kg) ~ log(Solid_P_mg_kg), data = subset(data_solids, 
##     Land_Coverage_Category == "Douglas-fir forest"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.12452 -0.24370  0.00926  0.39848  0.64306 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         14.2279     1.6388   8.682 1.61e-06 ***
## log(Solid_P_mg_kg)  -0.2846     0.2089  -1.363    0.198    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5105 on 12 degrees of freedom
## Multiple R-squared:  0.134,  Adjusted R-squared:  0.06181 
## F-statistic: 1.856 on 1 and 12 DF,  p-value: 0.1981
#sagebrush
model<-lm(log(Solid_C_mg_kg)~log(Solid_P_mg_kg), data = subset(data_solids, Land_Coverage_Category=="Sagebrush shrubland"))
summary(model)
## 
## Call:
## lm(formula = log(Solid_C_mg_kg) ~ log(Solid_P_mg_kg), data = subset(data_solids, 
##     Land_Coverage_Category == "Sagebrush shrubland"))
## 
## Residuals:
##        1        2        3        4        5 
## -0.18263  0.07222 -0.17107  0.19096  0.09052 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        14.76017    0.82811  17.824 0.000385 ***
## log(Solid_P_mg_kg) -0.29516    0.09526  -3.099 0.053357 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1936 on 3 degrees of freedom
## Multiple R-squared:  0.7619, Adjusted R-squared:  0.6826 
## F-statistic: 9.601 on 1 and 3 DF,  p-value: 0.05336
#P and N
ggplot(data_solids, aes(x=Solid_N_mg_kg, y=Solid_P_mg_kg, color = Land_Coverage_Category)) + geom_point() + geom_smooth(method="lm", se = FALSE) + theme_classic() + scale_color_manual(values = c("#346648", "#dbb40d"))
## `geom_smooth()` using formula = 'y ~ x'

#all feedstocks
model<-lm(log(Solid_N_mg_kg)~log(Solid_P_mg_kg), data = data_solids)
summary(model)
## 
## Call:
## lm(formula = log(Solid_N_mg_kg) ~ log(Solid_P_mg_kg), data = data_solids)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.84610 -0.30184  0.02675  0.34999  0.78060 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         9.08078    1.06949   8.491 1.61e-07 ***
## log(Solid_P_mg_kg) -0.05138    0.13241  -0.388    0.703    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4706 on 17 degrees of freedom
## Multiple R-squared:  0.008781,   Adjusted R-squared:  -0.04953 
## F-statistic: 0.1506 on 1 and 17 DF,  p-value: 0.7028
#Douglas fir only
model<-lm(log(Solid_N_mg_kg)~log(Solid_P_mg_kg), data = subset(data_solids, Land_Coverage_Category=="Douglas-fir forest"))
summary(model)
## 
## Call:
## lm(formula = log(Solid_N_mg_kg) ~ log(Solid_P_mg_kg), data = subset(data_solids, 
##     Land_Coverage_Category == "Douglas-fir forest"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.78941 -0.30421 -0.04207  0.37311  0.84806 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          9.4408     1.6029   5.890 7.37e-05 ***
## log(Solid_P_mg_kg)  -0.1090     0.2043  -0.534    0.603    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4993 on 12 degrees of freedom
## Multiple R-squared:  0.02318,    Adjusted R-squared:  -0.05822 
## F-statistic: 0.2848 on 1 and 12 DF,  p-value: 0.6033
#sagebrush
model<-lm(log(Solid_N_mg_kg)~log(Solid_P_mg_kg), data = subset(data_solids, Land_Coverage_Category=="Sagebrush shrubland"))
summary(model)
## 
## Call:
## lm(formula = log(Solid_N_mg_kg) ~ log(Solid_P_mg_kg), data = subset(data_solids, 
##     Land_Coverage_Category == "Sagebrush shrubland"))
## 
## Residuals:
##        1        2        3        4        5 
## -0.02507  0.31641  0.11150 -0.21428 -0.18856 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)   
## (Intercept)         10.8132     1.0893   9.927  0.00217 **
## log(Solid_P_mg_kg)  -0.2224     0.1253  -1.775  0.17400   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2547 on 3 degrees of freedom
## Multiple R-squared:  0.5122, Adjusted R-squared:  0.3496 
## F-statistic:  3.15 on 1 and 3 DF,  p-value: 0.174

####PCA & RDAs

#PCA of elemental concentration of the solid chars with burn severity and feedstock

#format dataframe using same code as correlation and then also remove NAs
data_solids_corr <- data_solids %>%
  rename_with(
    ~sub("Solid_(.*?)_mg_kg", "\\1", .),
    starts_with("Solid_")) %>%
  select(Burn_Severity, Land_Coverage_Category, P, C, N, Ca, Mg, Na, Fe, Al, S, K) %>%
  na.omit(data_solids)

pca <- prcomp(data_solids_corr[3:12], center=TRUE, scale=TRUE)

##Extract PCA information, including loadings, sample scores, and component percentages
eigval <- pca$sdev^2
eigvec <- data.frame(pca$rotation)
loadings1 <- eigvec$PC1*sqrt(eigval[1])
loadings2 <- eigvec$PC2*sqrt(eigval[2])
loadings <- cbind(loadings1, loadings2)

colnames(loadings)<-c('PC1','PC2')
rownames(loadings)<-colnames(data_solids_corr[3:12])
loadingstest<-as.data.frame(loadings)
pca_out <- as.data.frame(pca$x)
percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))


xscale <- 10
yscale <- 10
ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=data_solids_corr$Land_Coverage_Category,shape=data_solids_corr$Burn_Severity),size=3)+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)

#include pH, particulate concentrations from ICP, and filt0.7 concentrations from ICP, NPOC, TDN, MBD 

#fix THIS DOESN"T WORK RIGHT NOW
#format dataframe and then also remove NAs
#data_solids_corr <- data_solids %>%
#  rename_with(
#    ~sub("Leachate_(.*?)_mg_kg", "\\1", .),
#    starts_with("Leachate_")) %>%
#  select(Burn_Severity, Land_Coverage_Category, P, C, N, Ca, Mg, Na, Fe, Al, S, K) %>%
#  na.omit(data_solids)
#PCAs of XANES LCF

#FIX THIS DOENS'T RUN RIGHT NOW

#all phases together; open air; Po individual species

#format dataframe using same code as correlation and then also remove NAs
#XANES_corr <- XANES_clean %>%
#  select(Parent_ID, Burn_Severity, Land_Coverage_Category, FilterSize, XANES_Pi_Al, XANES_Pi_Ca, XANES_Pi_Fe, #XANES_Pi_K, XANES_Pi_Mg, XANES_Pi_Na, XANES_Po_Al, XANES_Po_Fe, XANES_Po_Ca, XANES_Po_Na) %>%
#  na.omit(XANES_clean)

#pca <- prcomp(XANES_corr[5:14], center=TRUE, scale=TRUE)

##Extract PCA information, including loadings, sample scores, and component percentages
#eigval <- pca$sdev^2
#eigvec <- data.frame(pca$rotation)
#loadings1 <- eigvec$PC1*sqrt(eigval[1])
#loadings2 <- eigvec$PC2*sqrt(eigval[2])
#loadings <- cbind(loadings1, loadings2)

#colnames(loadings)<-c('PC1','PC2')
#rownames(loadings)<-colnames(XANES_corr[5:14])
#loadingstest<-as.data.frame(loadings)
#pca_out <- as.data.frame(pca$x)
#percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
#percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))


#xscale <- 4
#yscale <- 4
#ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=XANES_corr$Land_Coverage_Category,shape=XANES_corr$Burn_Severity, size=XANES_corr$FilterSize))+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+scale_size_manual(values=c(3,6,9))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)+ggtitle("Open Air XANES - all phases")


#all phases together; open air; Po total instead of individual

#format dataframe using same code as correlation and then also remove NAs
#XANES_corr <- XANES_clean %>%
#  select(Parent_ID, Burn_Severity, Land_Coverage_Category, FilterSize, XANES_Pi_Al, XANES_Pi_Ca, XANES_Pi_Fe, XANES_Pi_K, XANES_Pi_Mg, XANES_Pi_Na, XANES_Po) %>%
#  na.omit(XANES_clean)

#pca <- prcomp(XANES_corr[5:11], center=TRUE, scale=TRUE)

##Extract PCA information, including loadings, sample scores, and component percentages
#eigval <- pca$sdev^2
#eigvec <- data.frame(pca$rotation)
#loadings1 <- eigvec$PC1*sqrt(eigval[1])
#loadings2 <- eigvec$PC2*sqrt(eigval[2])
#loadings <- cbind(loadings1, loadings2)

#colnames(loadings)<-c('PC1','PC2')
#rownames(loadings)<-colnames(XANES_corr[5:11])
#loadingstest<-as.data.frame(loadings)
#pca_out <- as.data.frame(pca$x)
#percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
#percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))


#xscale <- 4
#yscale <- 4
#ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=XANES_corr$Land_Coverage_Category,shape=XANES_corr$Burn_Severity, size=XANES_corr$FilterSize))+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+scale_size_manual(values=c(3,6,9))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)+ggtitle("Open Air XANES - all phases")



#solid only; open air; Po individual species

#format dataframe using same code as correlation and then also remove NAs; also removed columns that were all zeros for this phase
XANES_corr <- XANES_clean %>%
  select(Parent_ID, Burn_Severity, Land_Coverage_Category, FilterSize, XANES_Pi_Al, XANES_Pi_Ca, XANES_Pi_Fe, XANES_Pi_K, XANES_Pi_Mg, XANES_Pi_Na, XANES_Po_Ca, XANES_Po_Na) %>%
  na.omit(XANES_clean) %>%
  filter(FilterSize %in% c("solid"))

pca <- prcomp(XANES_corr[5:12], center=TRUE, scale=TRUE)

##Extract PCA information, including loadings, sample scores, and component percentages
eigval <- pca$sdev^2
eigvec <- data.frame(pca$rotation)
loadings1 <- eigvec$PC1*sqrt(eigval[1])
loadings2 <- eigvec$PC2*sqrt(eigval[2])
loadings <- cbind(loadings1, loadings2)

colnames(loadings)<-c('PC1','PC2')
rownames(loadings)<-colnames(XANES_corr[5:12])
loadingstest<-as.data.frame(loadings)
pca_out <- as.data.frame(pca$x)
percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))


xscale <- 4
yscale <- 4
ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=XANES_corr$Land_Coverage_Category,shape=XANES_corr$Burn_Severity),size=5)+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)+ggtitle("Open Air XANES - Solid")

#solid only; open air; Po grouped

#format dataframe using same code as correlation and then also remove NAs; also removed columns that were all zeros for this phase
XANES_corr <- XANES_clean %>%
  select(Parent_ID, Burn_Severity, Land_Coverage_Category, FilterSize, XANES_Pi_Al, XANES_Pi_Ca, XANES_Pi_Fe, XANES_Pi_K, XANES_Pi_Mg, XANES_Pi_Na, XANES_Po) %>%
  na.omit(XANES_clean) %>%
  filter(FilterSize %in% c("solid"))

pca <- prcomp(XANES_corr[5:11], center=TRUE, scale=TRUE)

##Extract PCA information, including loadings, sample scores, and component percentages
eigval <- pca$sdev^2
eigvec <- data.frame(pca$rotation)
loadings1 <- eigvec$PC1*sqrt(eigval[1])
loadings2 <- eigvec$PC2*sqrt(eigval[2])
loadings <- cbind(loadings1, loadings2)

colnames(loadings)<-c('PC1','PC2')
rownames(loadings)<-colnames(XANES_corr[5:11])
loadingstest<-as.data.frame(loadings)
pca_out <- as.data.frame(pca$x)
percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))


xscale <- 4
yscale <- 4
ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=XANES_corr$Land_Coverage_Category,shape=XANES_corr$Burn_Severity),size=5)+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)+ggtitle("Open Air XANES - Solid")

leachate <- data %>%
  select(1:3, 5, 10, 11, 49, 62, 70, 78, 50)
leachate_mean <- leachate %>%
  dplyr::group_by(Parent_ID) %>%
  dplyr::summarize(
    Leachate_P_mg_kgchar_unfilt_Mean = mean(Leachate_P_mg_kgchar_unfilt, na.rm = TRUE),
    Leachate_P_mg_kgchar_part_Mean = mean(Leachate_P_mg_kgchar_part, na.rm = TRUE),
    Leachate_P_mg_kgchar_filt0.7_Mean = mean(Leachate_P_mg_kgchar_filt0.7, na.rm = TRUE),
    MBD_P_mg_kgchar_Mean = mean(MBD_P_mg_kgchar, na.rm = TRUE),
    pH_Mean = mean(pH)
  )

solid <- data_solids %>%
  select(1,3, 5, 10, 11, 31:34, 88, 96:101, 111:120)

data_rda<- solid %>%
  full_join(leachate_mean) #%>%
## Joining with `by = join_by(Parent_ID)`
  #na.omit(solid) #remove NAs

data_rda <- data_rda %>%
  select(2,3,6:17, 27:30)

pca <- prcomp(data_rda[3:18], center=TRUE, scale=TRUE)

##Extract PCA information, including loadings, sample scores, and component percentages
eigval <- pca$sdev^2
eigvec <- data.frame(pca$rotation)
loadings1 <- eigvec$PC1*sqrt(eigval[1])
loadings2 <- eigvec$PC2*sqrt(eigval[2])
loadings <- cbind(loadings1, loadings2)

colnames(loadings)<-c('PC1','PC2')
rownames(loadings)<-colnames(data_rda[3:18])
loadingstest<-as.data.frame(loadings)
pca_out <- as.data.frame(pca$x)
percentage <- round(pca$sdev^2/sum(pca$sdev^2)*100,1)
percentage <- paste(colnames(pca_out), "-", paste(as.character(percentage), "%", sep=""))


xscale <- 6
yscale <- 6
ggplot()+geom_point(data=pca_out,mapping=aes(x=PC1, y=PC2, colour=data_rda$Land_Coverage_Category,shape=data_rda$Burn_Severity),size=5)+xlim(-xscale,xscale)+ylim(-yscale,yscale)+geom_abline(intercept=0, slope=0, linetype="dashed", size=0.8, colour="gray")+geom_vline(aes(xintercept=0), linetype="dashed", size=0.8, colour="gray")+theme_light()+labs(colour="Land_Coverage_Category",shape="Burn_Severity")+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),panel.background=element_blank(),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16),axis.text.y=element_text(size=16),axis.text.x=element_text(size=16),legend.title=element_text(size=12),legend.text=element_text(size=12),legend.position="bottom",legend.box="verticle",legend.margin=margin())+xlab(percentage[1])+ylab(percentage[2])+guides(color=guide_legend(nrow=2, byrow=TRUE))+scale_color_manual(values=c("#346648", "#dbb40d"))+scale_shape_manual(values=c(16,17,3,12))+geom_text_repel(data=loadingstest,mapping=aes(label=row.names(loadingstest), x=PC1*xscale, y=PC2*yscale),size=5, colour="black", alpha=0.5)

#setting up an RDA to look at how burn conditions (severity, duration, max temp) and chemistry of solid char (P, Ca, Mg, K, Fe, Al, Na concentration, XANES categories, and NMR categories) impact P mobilization (total P concentration in unfiltered leachate, aqueous phase (<0.7), and particulate phase (>0.7); soluble reactive P)


#this is example code of how to do an RDA
#reference: http://dmcglinn.github.io/quant_methods/lessons/multivariate_models.html
#other reference: https://popgen.nescent.org/2018-03-27_RDA_GEA.html
#citation in here: https://cran.r-project.org/web/packages/vegan/vignettes/FAQ-vegan.html
#reference about scaling: https://mb3is.megx.net/gustame/constrained-analyses/rda

#make dataframe with all leachate samples
leachate <- data %>%
  select(1:3, 5, 10, 11, 49, 62, 70, 78, 50)

solid <- data_solids %>%
  select(1, 31:34, 88, 96:101, 111:120)

data_rda<- leachate %>%
  full_join(solid) %>%
  filter(Solid_Na_g_kg != "") #remove NAs
## Joining with `by = join_by(Parent_ID)`
#create dataframe for sample variables
samples<-data_rda[ c(7:32)]
#create dataframe for burn condition variables
enviro<-data_rda[ c(4:6)]

#make dataframe with solid samples (leachate avg) to prevent using same number so many times
#alternative RDA: avergae leachate concentrations
leachate_mean <- leachate %>%
  dplyr::group_by(Parent_ID) %>%
  dplyr::summarize(
    Leachate_P_mg_kgchar_unfilt_Mean = mean(Leachate_P_mg_kgchar_unfilt, na.rm = TRUE),
    Leachate_P_mg_kgchar_part_Mean = mean(Leachate_P_mg_kgchar_part, na.rm = TRUE),
    Leachate_P_mg_kgchar_filt0.7_Mean = mean(Leachate_P_mg_kgchar_filt0.7, na.rm = TRUE),
    MBD_P_mg_kgchar_Mean = mean(MBD_P_mg_kgchar, na.rm = TRUE),
    pH_Mean = mean(pH)
  )

solid <- data_solids %>%
  select(1,3, 5, 10, 11, 31:34, 88, 96:101, 111:120)

data_rda<- solid %>%
  full_join(leachate_mean) %>%
  filter(Solid_Na_g_kg != "") #remove NAs
## Joining with `by = join_by(Parent_ID)`
#RDA: How do burn conditions (severity, duration, temp) influence chemistry of solids and leachates?
#create dataframe for sample variables
samples<-data_rda[ c(6:31)]
#create dataframe for burn condition variables
enviro<-data_rda[ c(3:5)]

#make sure both dataframes have equal rows
all.equal(rownames(samples), rownames(enviro))
## [1] TRUE
#RDA (this anaysis expects a linear response of each sample to the environmental variables)
rda_tree=rda(samples~ . , data=enviro)
rda_tree
## Call: rda(formula = samples ~ Burn_Severity + Burn_Duration +
## Char_Max_Temp, data = enviro)
## 
##                 Inertia Proportion Rank
## Total         1.975e+07  1.000e+00     
## Constrained   7.086e+06  3.588e-01    5
## Unconstrained 1.267e+07  6.412e-01   12
## Inertia is variance 
## 
## Eigenvalues for constrained axes:
##    RDA1    RDA2    RDA3    RDA4    RDA5 
## 7077630    5994    2089     345      87 
## 
## Eigenvalues for unconstrained axes:
##      PC1      PC2      PC3      PC4      PC5      PC6      PC7      PC8 
## 12637555    26075     1237      466      212      133       77       59 
##      PC9     PC10     PC11     PC12 
##       20        5        1        1
RsquareAdj(rda_tree)
## $r.squared
## [1] 0.358756
## 
## $adj.r.squared
## [1] 0.09157107
summary(eigenvals(rda_tree, model = "constrained"))
## Importance of components:
##                            RDA1      RDA2      RDA3      RDA4      RDA5
## Eigenvalue            7.078e+06 5.994e+03 2.089e+03 3.454e+02 8.651e+01
## Proportion Explained  9.988e-01 8.458e-04 2.948e-04 4.875e-05 1.221e-05
## Cumulative Proportion 9.988e-01 9.996e-01 9.999e-01 1.000e+00 1.000e+00
#see the factor loadings (correlations of each variable with each pca axis)
scores(rda_tree, choices=1:2, display="species", scaling=2)
##                                            RDA1          RDA2
## Ortho                              -0.258608414 -0.3836313766
## Mono                                0.174726143  0.3077998691
## Di                                  0.058838394  0.1025543877
## Pyro                                0.025383379 -0.0272017518
## XANES_Po                            0.281408979  0.5249502154
## XANES_Pi_Fe                         0.118469141  0.3483344631
## XANES_Pi_Na                         0.209494308 -0.2262644254
## XANES_Pi_Ca                        -0.248376554 -0.4292863568
## XANES_Pi_Al                        -0.011770369 -0.0002355077
## XANES_Pi_Mg                        -0.317177075 -0.1469297923
## XANES_Pi_K                         -0.014982768 -0.0554861412
## Solid_P_g_kg                       -0.080029673 -0.0212126139
## Solid_Ca_g_kg                      -0.365725846 -0.3100954065
## Solid_Fe_g_kg                      -0.050220016 -0.0196842734
## Solid_Al_g_kg                      -0.047826939 -0.0208178893
## Solid_K_g_kg                       -0.558503286 -0.0876518422
## Solid_Mg_g_kg                      -0.073643548 -0.0213093323
## Solid_Na_g_kg                      -0.009612414 -0.0046855763
## Solid_S_g_kg                       -0.021646259  0.0001316653
## Solid_C_g_kg                        1.091870534  2.1034905286
## Solid_N_g_kg                        0.019330517  0.0613465720
## Leachate_P_mg_kgchar_unfilt_Mean  -56.778548158  0.1556707251
## Leachate_P_mg_kgchar_part_Mean    -57.781948472 -0.0964623423
## Leachate_P_mg_kgchar_filt0.7_Mean   1.003400315  0.2521330673
## MBD_P_mg_kgchar_Mean                0.684190842  0.1321200187
## pH_Mean                            -0.050300766 -0.0401592675
## attr(,"const")
## [1] 135.3677
#plot it
plot(rda_tree, type='n', scaling=2)
orditorp(rda_tree, display='sp', cex=0.5, scaling=1, col='blue')
text(rda_tree, display='cn', col='red')

#check correlation
pairs.panels(enviro, scale=T)

#can visualize RDA importance
screeplot(rda_tree)

#check significance of all PCAs
signif.full <- anova.cca(rda_tree, parallel=getOption("mc.cores")) # default is permutation=999
signif.full #model is significant
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Burn_Severity + Burn_Duration + Char_Max_Temp, data = enviro)
##          Df Variance      F Pr(>F)
## Model     5  7086145 1.3427  0.259
## Residual 12 12665842
#check significance of individual PCAs
signif.axis <- anova.cca(rda_tree, by="axis", parallel=getOption("mc.cores"))
signif.axis 
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Burn_Severity + Burn_Duration + Char_Max_Temp, data = enviro)
##          Df Variance      F Pr(>F)
## RDA1      1  7077630 6.7056  0.272
## RDA2      1     5994 0.0057  1.000
## RDA3      1     2089 0.0020  1.000
## RDA4      1      345 0.0003  1.000
## RDA5      1       87 0.0001  1.000
## Residual 12 12665842
#check VIF
vif.cca(rda_tree)
##      Burn_SeverityLow Burn_SeverityModerate Burn_SeverityUnburned 
##              3.437847              2.249385              8.802287 
##         Burn_Duration         Char_Max_Temp 
##              2.638543              8.735625
#compare the model to chance with permutation test
anova(rda_tree)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Burn_Severity + Burn_Duration + Char_Max_Temp, data = enviro)
##          Df Variance      F Pr(>F)
## Model     5  7086145 1.3427   0.28
## Residual 12 12665842
#plot again
plot(rda_tree, scaling=2) #this just doesn't work; not significant, VIF too high

#RDA: How does solid P chemistry (NMR, XANES) influence elemental concentrations in leachate (unfilt, filt, part, molybdate)?; color by burn severity, shape by feedstock; uses leachate means; VIF is too high
data_rda<- solid %>%
  full_join(leachate_mean) %>%
  filter(Solid_Na_g_kg != "") #remove NAs
## Joining with `by = join_by(Parent_ID)`
#create dataframe for sample variables
samples<-data_rda[ c(27:30)]
#create dataframe for burn condition variables
enviro<-data_rda[ c(6:16)]

#make sure both dataframes have equal rows
all.equal(rownames(samples), rownames(enviro))
## [1] TRUE
#RDA (this anaysis expects a linear response of each sample to the environmental variables)
rda_tree=rda(samples~ . , data=enviro)
rda_tree
## Call: rda(formula = samples ~ Ortho + Mono + Di + Pyro + XANES_Po +
## XANES_Pi_Fe + XANES_Pi_Na + XANES_Pi_Ca + XANES_Pi_Al + XANES_Pi_Mg +
## XANES_Pi_K, data = enviro)
## 
##                 Inertia Proportion Rank
## Total         1.974e+07  1.000e+00     
## Constrained   1.651e+07  8.366e-01    3
## Unconstrained 3.225e+06  1.634e-01    3
## Inertia is variance 
## 
## Eigenvalues for constrained axes:
##     RDA1     RDA2     RDA3 
## 16497551    16809      359 
## 
## Eigenvalues for unconstrained axes:
##     PC1     PC2     PC3 
## 3213434   12019      38
RsquareAdj(rda_tree)
## $r.squared
## [1] 0.836603
## 
## $adj.r.squared
## [1] 0.5370419
summary(eigenvals(rda_tree, model = "constrained"))
## Importance of components:
##                           RDA1      RDA2      RDA3
## Eigenvalue            1.65e+07 1.681e+04 3.591e+02
## Proportion Explained  9.99e-01 1.018e-03 2.174e-05
## Cumulative Proportion 9.99e-01 1.000e+00 1.000e+00
#see the factor loadings (correlations of each variable with each pca axis)
scores(rda_tree, choices=1:2, display="species", scaling=2)
##                                         RDA1      RDA2
## Leachate_P_mg_kgchar_unfilt_Mean  87.3227700 -1.432921
## Leachate_P_mg_kgchar_part_Mean    87.6604394  1.413790
## Leachate_P_mg_kgchar_filt0.7_Mean -0.3376693 -2.846711
## MBD_P_mg_kgchar_Mean              -0.1250006 -1.855479
## attr(,"const")
## [1] 135.3475
#plot it
plot(rda_tree, type='n', scaling=2)
orditorp(rda_tree, display='sp', cex=0.5, scaling=1, col='blue')
text(rda_tree, display='cn', col='red')

#check correlation
pairs.panels(enviro, scale=T)

#can visualize RDA importance
screeplot(rda_tree)

#check significance of all PCAs
signif.full <- anova.cca(rda_tree, parallel=getOption("mc.cores")) # default is permutation=999
signif.full #model is significant
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Ortho + Mono + Di + Pyro + XANES_Po + XANES_Pi_Fe + XANES_Pi_Na + XANES_Pi_Ca + XANES_Pi_Al + XANES_Pi_Mg + XANES_Pi_K, data = enviro)
##          Df Variance      F Pr(>F)
## Model    11 16514719 2.7928  0.133
## Residual  6  3225491
#check significance of individual PCAs
signif.axis <- anova.cca(rda_tree, by="axis", parallel=getOption("mc.cores"))
signif.axis #only RDA1 and RDA 2 are sig
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Ortho + Mono + Di + Pyro + XANES_Po + XANES_Pi_Fe + XANES_Pi_Na + XANES_Pi_Ca + XANES_Pi_Al + XANES_Pi_Mg + XANES_Pi_K, data = enviro)
##          Df Variance       F Pr(>F)
## RDA1      1 16497551 71.6064  0.121
## RDA2      1    16809  0.0730  1.000
## RDA3      1      359  0.0016  1.000
## Residual 14  3225491
#check VIF
vif.cca(rda_tree)
##        Ortho         Mono           Di         Pyro     XANES_Po  XANES_Pi_Fe 
## 430897.99439 244137.82681  27453.58296  31271.68532   2066.14942   1017.06394 
##  XANES_Pi_Na  XANES_Pi_Ca  XANES_Pi_Al  XANES_Pi_Mg   XANES_Pi_K 
##   2069.52735   1022.25453     83.70275   1425.22502    481.65760
#compare the model to chance with permutation test
anova(rda_tree)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Ortho + Mono + Di + Pyro + XANES_Po + XANES_Pi_Fe + XANES_Pi_Na + XANES_Pi_Ca + XANES_Pi_Al + XANES_Pi_Mg + XANES_Pi_K, data = enviro)
##          Df Variance      F Pr(>F)
## Model    11 16514719 2.7928  0.135
## Residual  6  3225491
#plot again
plot(rda_tree, scaling=2) #this just doesn't work; not significant, VIF too high

#How does feedstock and burn severity influence solid P chemistry (P conc, NMR, XANES)?
data_rda<- solid %>%
  full_join(leachate_mean) 
## Joining with `by = join_by(Parent_ID)`
#create dataframe for sample variables
samples<-data_rda[ c(6:17)]
#create dataframe for burn condition variables
enviro<-data_rda[ c(2:3)]

#make sure both dataframes have equal rows
all.equal(rownames(samples), rownames(enviro))
## [1] TRUE
#RDA (this anaysis expects a linear response of each sample to the environmental variables)
rda_tree=rda(samples~ . , data=enviro)
rda_tree
## Call: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity,
## data = enviro)
## 
##                 Inertia Proportion Rank
## Total         2544.1397     1.0000     
## Constrained   1767.1852     0.6946    4
## Unconstrained  776.9546     0.3054   12
## Inertia is variance 
## 
## Eigenvalues for constrained axes:
##   RDA1   RDA2   RDA3   RDA4 
## 1300.3  399.3   37.0   30.6 
## 
## Eigenvalues for unconstrained axes:
##   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11  PC12 
## 420.4 147.8  83.8  71.1  31.6  16.2   2.9   2.3   0.7   0.3   0.0   0.0
RsquareAdj(rda_tree)
## $r.squared
## [1] 0.6946101
## 
## $adj.r.squared
## [1] 0.6073559
summary(eigenvals(rda_tree, model = "constrained"))
## Importance of components:
##                            RDA1     RDA2     RDA3     RDA4
## Eigenvalue            1300.2585 399.2663 37.02949 30.63083
## Proportion Explained     0.7358   0.2259  0.02095  0.01733
## Cumulative Proportion    0.7358   0.9617  0.98267  1.00000
#see the factor loadings (correlations of each variable with each pca axis)
scores(rda_tree, choices=1:2, display="species", scaling=2)
##                     RDA1       RDA2
## Ortho         4.12832420  1.1490159
## Mono         -3.24044129 -0.2501559
## Di           -1.20235534  0.2376141
## Pyro          0.31901678 -1.1398120
## XANES_Po     -5.78118619  0.9835438
## XANES_Pi_Fe  -3.73794387  1.3954752
## XANES_Pi_Na   0.88496193 -4.8047389
## XANES_Pi_Ca   4.49340039  0.7683397
## XANES_Pi_Al   0.08413426 -0.4895688
## XANES_Pi_Mg   3.35883295  1.8850168
## XANES_Pi_K    0.50231586  0.2399660
## Solid_P_g_kg  0.60750473  0.6357070
## attr(,"const")
## [1] 14.62862
#plot it
plot(rda_tree, type='n', scaling=2)
orditorp(rda_tree, display='sp', cex=0.5, scaling=1, col='blue')
text(rda_tree, display='cn', col='red')

#check correlation
pairs.panels(enviro, scale=T)

#can visualize RDA importance
screeplot(rda_tree)

#check significance of all PCAs
signif.full <- anova.cca(rda_tree, parallel=getOption("mc.cores")) # default is permutation=999
signif.full #model is significant
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
##          Df Variance      F Pr(>F)    
## Model     4  1767.19 7.9608  0.001 ***
## Residual 14   776.95                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#check significance of individual PCAs
signif.axis <- anova.cca(rda_tree, by="axis", parallel=getOption("mc.cores"))
signif.axis #only RDA1 and RDA 2 are sig
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
##          Df Variance       F Pr(>F)    
## RDA1      1  1300.26 23.4295  0.001 ***
## RDA2      1   399.27  7.1944  0.014 *  
## RDA3      1    37.03  0.6672  0.903    
## RDA4      1    30.63  0.5519  0.707    
## Residual 14   776.95                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#check VIF
vif.cca(rda_tree)
## Land_Coverage_CategorySagebrush shrubland 
##                                  1.118041 
##                          Burn_SeverityLow 
##                                  1.846364 
##                     Burn_SeverityModerate 
##                                  1.836781 
##                     Burn_SeverityUnburned 
##                                  1.558868
#compare the model to chance with permutation test
anova(rda_tree)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
##          Df Variance      F Pr(>F)    
## Model     4  1767.19 7.9608  0.001 ***
## Residual 14   776.95                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#plot again
plot(rda_tree, scaling=2) #this is significant and VIF is good; make this one pretty

#How does feedstock and burn severity influence solid (P concentration, NMR, XANES) and leachate (concentration of different phases and molybdate) chemistry 
data_rda<- solid %>%
  full_join(leachate_mean) 
## Joining with `by = join_by(Parent_ID)`
#create dataframe for sample variables
samples<-data_rda[ c(6:17, 27:30)]
#create dataframe for burn condition variables
enviro<-data_rda[ c(2:3)]

#make sure both dataframes have equal rows
all.equal(rownames(samples), rownames(enviro))
## [1] TRUE
#RDA (this anaysis expects a linear response of each sample to the environmental variables)
rda_tree=rda(samples~ . , data=enviro)
rda_tree
## Call: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity,
## data = enviro)
## 
##                 Inertia Proportion Rank
## Total         1.882e+07  1.000e+00     
## Constrained   7.940e+06  4.219e-01    4
## Unconstrained 1.088e+07  5.781e-01   12
## Inertia is variance 
## 
## Eigenvalues for constrained axes:
##    RDA1    RDA2    RDA3    RDA4 
## 7926038   12958     984      81 
## 
## Eigenvalues for unconstrained axes:
##      PC1      PC2      PC3      PC4      PC5      PC6      PC7      PC8 
## 10869073    11608      426      214       92       75       48       17 
##      PC9     PC10     PC11     PC12 
##        8        2        1        0
RsquareAdj(rda_tree)
## $r.squared
## [1] 0.4218584
## 
## $adj.r.squared
## [1] 0.2566751
summary(eigenvals(rda_tree, model = "constrained"))
## Importance of components:
##                            RDA1      RDA2      RDA3      RDA4
## Eigenvalue            7.926e+06 1.296e+04 9.840e+02 8.051e+01
## Proportion Explained  9.982e-01 1.632e-03 1.239e-04 1.014e-05
## Cumulative Proportion 9.982e-01 9.999e-01 1.000e+00 1.000e+00
#see the factor loadings (correlations of each variable with each pca axis)
scores(rda_tree, choices=1:2, display="species", scaling=2)
##                                          RDA1        RDA2
## Ortho                              0.31375232 -0.18493684
## Mono                              -0.19173524  0.13469572
## Di                                -0.04856659  0.03015374
## Pyro                              -0.07313601  0.02074476
## XANES_Po                          -0.27017217  0.03851785
## XANES_Pi_Fe                       -0.14214762  0.01074035
## XANES_Pi_Na                       -0.31957053  0.07374803
## XANES_Pi_Ca                        0.29969687 -0.14665261
## XANES_Pi_Al                       -0.03611682 -0.01006229
## XANES_Pi_Mg                        0.39900322 -0.05409895
## XANES_Pi_K                         0.05464259  0.09400254
## Solid_P_g_kg                       0.10505008  0.01294289
## Leachate_P_mg_kgchar_unfilt_Mean  62.56938415  1.27600954
## Leachate_P_mg_kgchar_part_Mean    61.92711377 -1.32452727
## Leachate_P_mg_kgchar_filt0.7_Mean  0.64227038  2.60053681
## MBD_P_mg_kgchar_Mean               0.44612320  1.55961875
## attr(,"const")
## [1] 135.6696
#plot it
plot(rda_tree, type='n', scaling=2)
orditorp(rda_tree, display='sp', cex=0.5, scaling=1, col='blue')
text(rda_tree, display='cn', col='red')

#check correlation
pairs.panels(enviro, scale=T)

#can visualize RDA importance
screeplot(rda_tree)

#check significance of all PCAs
signif.full <- anova.cca(rda_tree, parallel=getOption("mc.cores")) # default is permutation=999
signif.full #model is significant
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
##          Df Variance      F Pr(>F)  
## Model     4  7940061 2.5539  0.059 .
## Residual 14 10881564                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#check significance of individual PCAs
signif.axis <- anova.cca(rda_tree, by="axis", parallel=getOption("mc.cores"))
signif.axis #only RDA1 and RDA 2 are sig
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
##          Df Variance       F Pr(>F)  
## RDA1      1  7926038 10.1975  0.058 .
## RDA2      1    12958  0.0167  0.999  
## RDA3      1      984  0.0013  1.000  
## RDA4      1       81  0.0001  0.999  
## Residual 14 10881564                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#check VIF
vif.cca(rda_tree)
## Land_Coverage_CategorySagebrush shrubland 
##                                  1.118041 
##                          Burn_SeverityLow 
##                                  1.846364 
##                     Burn_SeverityModerate 
##                                  1.836781 
##                     Burn_SeverityUnburned 
##                                  1.558868
#compare the model to chance with permutation test
anova(rda_tree)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = samples ~ Land_Coverage_Category + Burn_Severity, data = enviro)
##          Df Variance      F Pr(>F)  
## Model     4  7940061 2.5539  0.057 .
## Residual 14 10881564                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#plot again
plot(rda_tree, scaling=2) #this is isn't as good with the leachates included

#make a pretty plot
#use
#USE THIS ONE AS AN EXAMPLE!! (has horizons as difference shapes but need to figure out legend)
#levels(enviro$Site) <- c("GR2","GR3","GR5","BAR","SJER","SPR", "PR", "SH")
#eco <- enviro$Site
#bgg <- c("#000080", "#8B0000","#0E5E0B","#606060","#5B94BA","#DF1111", "#069710",  "#C9D4D1") # 8 colors for our sites
#shapes <- c(18, 16, 15,17) 
#sh <- as.factor(enviro$Horizon)
#plot(rda_tree, type="n", scaling=2)
#points(rda_tree, display="species", pch=20, cex=0.7, col="gray32", scaling=2)           #the chemical species from NMR and XANES
#points(rda_tree, display="sites", pch=shapes[sh], cex=1.3, scaling=2, col=bgg[eco]) # the soil samples
#text(rda_tree, scaling=2, display="bp", col="#0868ac", cex=1)                           # the predictors
#legend("bottomright", legend=levels(eco), bty="n", pch=21, cex=1, col=bg)
#legend("bottomright", legend=c("A", "B", "C", "GR2", "GR3", "GR5", "BAR", "SJER", "SPR", "PR", "SH"), col=c(rep("black",3), "#000080", "#8B0000","#0E5E0B","#606060","#5B94BA","#DF1111", "#069710",  "#C9D4D1"), pch=c(16, 15, 17, 16, 16, 16, 16, 16, 16, 16, 16))


#surface horizons only
#Master_CZO_WM_Soil_5.8.2020 <- read.csv("~/Documents/Research Notes/Dissertation General/Master_CZO_WM_Soil_5.8.2020.csv")
#Master_CZO_WM_Soil_5.8.2020$Climo<-factor(Master_CZO_WM_Soil_5.8.2020$Climo, levels=c("WM", "CZO"))
#xas.percent<-subset(Master_CZO_WM_Soil_5.8.2020, Sample.Type=="Soil")
#xas.percent<-xas.percent[ , c("Horizon.Category.C",  "Climo", "Site", "EEMT", "AI", "MAT", "MAP","pH.CaCl2", "Clay.Percent", "Fe.g.kgsoil.PMnorm", "Ca.g.kgsoil.PMnorm", "Al.g.kgsoil.PMnorm", "C.g.kgsoil", "P.g.kgsoil.PMnorm", "N.g.kgsoil", "Ca.Pi.Percent", "Al.Pi.Percent", "Fe.Pi.Percent", "XANES.Total.Monoester.Percent", "XANES.Total.Diester.Percent")]
#xas.percent<-subset(xas.percent, Ca.Pi.Percent!="NA")
#colnames(xas.percent)[1:20]<-c("Horizon",  "Gradient", "Site","EEMT", "Humidity", "MAT", "MAP","pH", "Clay", "Fe", "Ca", "Al", "C", "P", "N", "Ca-Pi", "Al-Pi", "Fe-Pi", "Monoester", "Diester")
#xas.percent<-subset(xas.percent, Horizon=="A")

####stacked spectra and pie charts

#example stacked spectra
NMR_spectra_df_subset <-NMR_spectra %>%
  select("Chemical_Shift_parts_per_million", "BSLE_0013-solid", "BSLE_0007-solid","BSLE_0002-solid", "BSLE_0050-solid") %>%
  dplyr::rename('solid_13' = 'BSLE_0013-solid',
                'solid_7' = 'BSLE_0007-solid',
                'solid_2' = 'BSLE_0002-solid',
                'solid_50' = 'BSLE_0050-solid') %>%
  mutate_all(~replace(., . == 0, NA)) %>%
  mutate(solid_13 = solid_13 + 1,
         solid_7 = solid_7 + 75,
         solid_2 = solid_2 + 150,
         solid_50 = solid_50 + 225) 
         
NMR_spectra_df_subset<-melt(NMR_spectra_df_subset, id=c("Chemical_Shift_parts_per_million"))

NMR_spectra_df_subset<-subset(NMR_spectra_df_subset, value!="")

NMR_stacked_spectra_df <- ggplot(NMR_spectra_df_subset, aes(x=Chemical_Shift_parts_per_million, y = value, color = as.factor(variable))) + geom_path()+theme_classic() + ylim(0,300) + scale_x_reverse() +xlim(10,-10) +scale_color_manual(values=c('Black','Black', 'Black', 'Black')) + theme(axis.ticks.y = element_blank(),  axis.text.y = element_blank(), legend.position = "none")  + ylab("") + xlab("Chemical Shift (ppm)")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
#cowplot::save_plot("../figures/NMR_stacked_spectra_df.pdf", NMR_stacked_spectra_df, ncol = 1, nrow = 1, base_aspect_ratio= 1:1.7, dpi=72)

pdf(file = "../figures/NMR_stacked_spectra_df.pdf", width = 3.5, height = 7) 
print(NMR_stacked_spectra_df)
## Warning: Removed 34204 rows containing missing values (`geom_path()`).
dev.off()
## quartz_off_screen 
##                 2
NMR_spectra_sb_subset <-NMR_spectra %>%
  select("Chemical_Shift_parts_per_million", "BSLE_0014-solid", "BSLE_0011-solid","BSLE_0072-solid") %>%
  dplyr::rename('solid_14' = 'BSLE_0014-solid',
                'solid_11' = 'BSLE_0011-solid',
                'solid_72' = 'BSLE_0072-solid') %>%
  mutate_all(~replace(., . == 0, NA)) %>%
  mutate(solid_14 = solid_14 + 1,
         solid_11 = solid_11 + 75,
         solid_72 = solid_72 + 150) 
         
NMR_spectra_sb_subset<-melt(NMR_spectra_sb_subset, id=c("Chemical_Shift_parts_per_million"))

NMR_spectra_sb_subset<-subset(NMR_spectra_sb_subset, value!="")

NMR_stacked_spectra_sb <- ggplot(NMR_spectra_sb_subset, aes(x=Chemical_Shift_parts_per_million, y = value, color = as.factor(variable))) + geom_path()+theme_classic() + ylim(0,300) + scale_x_reverse() +xlim(10,-10) +scale_color_manual(values=c('Black','Black', 'Black', 'Black')) + theme(axis.ticks.y = element_blank(),  axis.text.y = element_blank(), legend.position = "none")  + ylab("") + xlab("Chemical Shift (ppm)")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
#cowplot::save_plot("../figures/NMR_stacked_spectra_sb.pdf", NMR_stacked_spectra_sb, ncol = 1, nrow = 1, base_aspect_ratio= 1:1.7, dpi=72)

pdf(file = "../figures/NMR_stacked_spectra_sb.pdf", width = 3.5, height = 7) 
print(NMR_stacked_spectra_sb)
## Warning: Removed 25652 rows containing missing values (`geom_path()`).
dev.off()
## quartz_off_screen 
##                 2
#pie charts
NMR <- summary_solids_conc %>%
  select("Land_Coverage_Category", "Burn_Severity","Ortho_Mean", "Mono_Mean", "Di_Mean", "Pyro_Mean") %>%
  dplyr::rename(Ortho = Ortho_Mean,
                Mono = Mono_Mean,
                Di = Di_Mean,
                Pyro = Pyro_Mean)


NMR_melt<-melt(NMR, id=c("Land_Coverage_Category", "Burn_Severity"))

NMR_melt$variable=factor(NMR_melt$variable, levels= c("Ortho", "Pyro", "Mono", "Di"))

DF.raw<-subset(NMR_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Unburned")
nmr.legend<-ggplot(DF.raw, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
    theme(legend.key.size = unit(1, 'cm'), #change legend key size
          legend.key.height = unit(1, 'cm'),  legend.title = element_text(size=25),
          legend.key.width = unit(1, 'cm'), #change legend key width
          legend.text = element_text(size=20))+labs(fill="")
pdf(file = "../figures/nmr.legend.pdf", width = 7, height = 3.5) 
print(nmr.legend)
dev.off()
## quartz_off_screen 
##                 2
DF.raw.nmr<-ggplot(DF.raw, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/DF.raw.nmr.pdf", width = 1, height = 1) 
print(DF.raw.nmr)
dev.off()
## quartz_off_screen 
##                 2
DF.low<-subset(NMR_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Low")
DF.low.nmr<-ggplot(DF.low, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/DF.low.nmr.pdf", width = 1, height = 1) 
print(DF.low.nmr)
dev.off()
## quartz_off_screen 
##                 2
DF.mod<-subset(NMR_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Moderate")
DF.mod.nmr<-ggplot(DF.mod, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/DF.mod.nmr.pdf", width = 1, height = 1) 
print(DF.mod.nmr)
dev.off()
## quartz_off_screen 
##                 2
DF.high<-subset(NMR_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="High")
DF.high.nmr<-ggplot(DF.high, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/DF.high.nmr.pdf", width = 1, height = 1) 
print(DF.high.nmr)
dev.off()
## quartz_off_screen 
##                 2
SB.raw<-subset(NMR_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Unburned")
SB.raw.nmr<-ggplot(SB.raw, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/SB.raw.nmr.pdf", width = 1, height = 1) 
print(SB.raw.nmr)
dev.off()
## quartz_off_screen 
##                 2
SB.low<-subset(NMR_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Low")
SB.low.nmr<-ggplot(SB.low, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/SB.low.nmr.pdf", width = 1, height = 1) 
print(SB.low.nmr)
dev.off()
## quartz_off_screen 
##                 2
SB.mod<-subset(NMR_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Moderate")
SB.mod.nmr<-ggplot(SB.mod, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#780000", "#c1121f", "#003049", "#669BBC"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/SB.mod.nmr.pdf", width = 1, height = 1) 
print(SB.mod.nmr)
dev.off()
## quartz_off_screen 
##                 2
#example stacked spectra
XANES_spectra_df_subset <-XANES_spectra %>%
  select("eV", "BSLE_0013-solid_P-XANES.xmu.nor", "BSLE_0007-solid_P-XANES.xmu.nor","BSLE_0002-solid_P-XANES.xmu.nor", "BSLE_0050-solid_P-XANES.xmu.nor") %>%
  dplyr::rename('solid_13' = 'BSLE_0013-solid_P-XANES.xmu.nor',
                'solid_7' = 'BSLE_0007-solid_P-XANES.xmu.nor',
                'solid_2' = 'BSLE_0002-solid_P-XANES.xmu.nor',
                'solid_50' = 'BSLE_0050-solid_P-XANES.xmu.nor') %>%
  mutate_all(~replace(., . == 0, NA)) %>%
  mutate(solid_13 = solid_13 + 1,
         solid_7 = solid_7 + 8,
         solid_2 = solid_2 + 16,
         solid_50 = solid_50 + 24) 
         
XANES_spectra_df_subset<-melt(XANES_spectra_df_subset, id=c("eV"))

XANES_spectra_df_subset<-subset(XANES_spectra_df_subset, value!="")

XANES_stacked_spectra_df <- ggplot(XANES_spectra_df_subset, aes(x=eV, y = value, color = as.factor(variable))) + geom_path()+theme_classic() + scale_color_manual(values=c('Black','Black', 'Black', 'Black')) + theme(axis.ticks.y = element_blank(),  axis.text.y = element_blank(), legend.position = "none")  + ylab("") + xlab("Energy (eV)") + xlim(2140, 2200) + ylab ("Normalized Absorption (arb. units)")

pdf(file = "../figures/XANES_stacked_spectra_df.pdf", width = 3, height = 7) 
print(XANES_stacked_spectra_df)
## Warning: Removed 188 rows containing missing values (`geom_path()`).
dev.off()
## quartz_off_screen 
##                 2
XANES_spectra_sb_subset <-XANES_spectra %>%
  select("eV", "BSLE_0014-solid_P-XANES.xmu.nor", "BSLE_0011-solid_P-XANES.xmu.nor","BSLE_0072-solid_P-XANES.xmu.nor") %>%
  dplyr::rename('solid_14' = 'BSLE_0014-solid_P-XANES.xmu.nor',
                'solid_11' = 'BSLE_0011-solid_P-XANES.xmu.nor',
                'solid_72' = 'BSLE_0072-solid_P-XANES.xmu.nor') %>%
  mutate_all(~replace(., . == 0, NA)) %>%
  mutate(blank = solid_72) %>%
  mutate(solid_14 = solid_14 + 1,
         solid_11 = solid_11 + 8,
         solid_72 = solid_72 + 16,
         blank = blank + 24) 
         
XANES_spectra_sb_subset<-melt(XANES_spectra_sb_subset, id=c("eV"))

XANES_spectra_sb_subset<-subset(XANES_spectra_sb_subset, value!="")

XANES_stacked_spectra_sb <- ggplot(XANES_spectra_sb_subset, aes(x=eV, y = value, color = as.factor(variable))) + geom_path()+theme_classic() + scale_color_manual(values=c('Black','Black', 'Black', 'White')) + theme(axis.ticks.y = element_blank(),  axis.text.y = element_blank(), legend.position = "none")  + ylab("") + xlab("Energy (eV)") + xlim(2140, 2200) + ylab ("Normalized Absorption (arb. units)")

#cowplot::save_plot("../figures/NMR_stacked_spectra_sb.pdf", NMR_stacked_spectra_sb, ncol = 1, nrow = 1, base_aspect_ratio= 1:1.7, dpi=72)

pdf(file = "../figures/XANES_stacked_spectra_sb.pdf", width = 3, height = 7) 
print(XANES_stacked_spectra_sb)
## Warning: Removed 188 rows containing missing values (`geom_path()`).
dev.off()
## quartz_off_screen 
##                 2
#pie charts with Po grouped and inorganic separate metals
XANES <- summary_solids_conc %>%
  select("Land_Coverage_Category", "Burn_Severity","XANES_Pi_K_Mean", "XANES_Pi_Mg_Mean", "XANES_Pi_Al_Mean", "XANES_Pi_Ca_Mean", "XANES_Pi_Na_Mean", "XANES_Pi_Fe_Mean", "XANES_Po_Mean") %>%
  dplyr::rename(K = XANES_Pi_K_Mean,
                Mg = XANES_Pi_Mg_Mean,
                Al = XANES_Pi_Al_Mean,
                Ca = XANES_Pi_Ca_Mean,
                Na = XANES_Pi_Na_Mean,
                Fe = XANES_Pi_Fe_Mean,
                Po = XANES_Po_Mean)

XANES_melt<-melt(XANES, id=c("Land_Coverage_Category", "Burn_Severity"))

#XANES_melt$variable=factor(XANES_melt$variable, levels= c("Ortho", "Pyro", "Mono", "Di"))

DF.raw<-subset(XANES_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Unburned")
xanes.legend<-ggplot(DF.raw, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
    theme(legend.key.size = unit(1, 'cm'), #change legend key size
          legend.key.height = unit(1, 'cm'),  legend.title = element_text(size=25),
          legend.key.width = unit(1, 'cm'), #change legend key width
          legend.text = element_text(size=20))+labs(fill="")
pdf(file = "../figures/xanes.legend.pdf", width = 7, height = 3.5) 
print(xanes.legend)
dev.off()
## quartz_off_screen 
##                 2
DF.raw.xanes<-ggplot(DF.raw, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/DF.raw.xanes.pdf", width = 1, height = 1) 
print(DF.raw.xanes)
dev.off()
## quartz_off_screen 
##                 2
DF.low<-subset(XANES_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Low")
DF.low.xanes<-ggplot(DF.low, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/DF.low.xanes.pdf", width = 1, height = 1) 
print(DF.low.xanes)
dev.off()
## quartz_off_screen 
##                 2
DF.mod<-subset(XANES_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="Moderate")
DF.mod.xanes<-ggplot(DF.mod, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/DF.mod.xanes.pdf", width = 1, height = 1) 
print(DF.mod.xanes)
dev.off()
## quartz_off_screen 
##                 2
DF.high<-subset(XANES_melt, Land_Coverage_Category=="Douglas-fir forest" & Burn_Severity=="High")
DF.high.xanes<-ggplot(DF.high, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/DF.high.xanes.pdf", width = 1, height = 1) 
print(DF.high.xanes)
dev.off()
## quartz_off_screen 
##                 2
SB.raw<-subset(XANES_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Unburned")
SB.raw.xanes<-ggplot(SB.raw, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/SB.raw.xanes.pdf", width = 1, height = 1) 
print(SB.raw.xanes)
dev.off()
## quartz_off_screen 
##                 2
SB.low<-subset(XANES_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Low")
SB.low.xanes<-ggplot(SB.low, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/SB.low.xanes.pdf", width = 1, height = 1) 
print(SB.low.xanes)
dev.off()
## quartz_off_screen 
##                 2
SB.mod<-subset(XANES_melt, Land_Coverage_Category=="Sagebrush shrubland" & Burn_Severity=="Moderate")
SB.mod.xanes<-ggplot(SB.mod, aes(x = "", y = value, fill = variable)) +
    geom_col(color = "black") +
    coord_polar(theta = "y") + scale_fill_manual(values = c("#424a26", "#9B9B7A", "#D9AE94", "#F1DCA7", "#FFCB69", "#D08C60", "#997B66"))+theme_void()+
    theme(legend.position = "none")
pdf(file = "../figures/SB.mod.xanes.pdf", width = 1, height = 1) 
print(SB.mod.xanes)
dev.off()
## quartz_off_screen 
##                 2

#SEM

#reference: https://stats.oarc.ucla.edu/r/seminars/rsem/ #followed Model 4A

#make dataframe with solid samples (leachate avg) to prevent using same number so many times
leachate_mean <- leachate %>%
  dplyr::group_by(Parent_ID) %>%
  dplyr::summarize(
    Leachate_P_mg_kgchar_unfilt_Mean = mean(Leachate_P_mg_kgchar_unfilt, na.rm = TRUE),
    Leachate_P_mg_kgchar_part_Mean = mean(Leachate_P_mg_kgchar_part, na.rm = TRUE),
    Leachate_P_mg_kgchar_filt0.7_Mean = mean(Leachate_P_mg_kgchar_filt0.7, na.rm = TRUE),
    MBD_P_mg_kgchar_Mean = mean(MBD_P_mg_kgchar, na.rm = TRUE),
    pH_Mean = mean(pH)
  )

solid <- data_solids %>%
  select(1,3, 5, 10, 11, 31:34, 88, 96:101, 111:120)

data_rda<- solid %>%
  full_join(leachate_mean) %>%
  filter(Solid_Na_g_kg != "") #remove NAs
## Joining with `by = join_by(Parent_ID)`
model <- '
  # Measurement Model
    f1 =~ Land_Coverage_Category + Char_Max_Temp + Burn_Duration
    f2 =~ Burn_Severity

  # Structural Model
    f1 ~ a*f2
    y1 ~ b*f1

  # Residual Variances
    Land_Coverage_Category ~~ Land_Coverage_Category
    Char_Max_Temp ~~ Char_Max_Temp
    Burn_Duration ~~ Burn_Duration
    Burn_Severity ~~ Burn_Severity
'

Leachate_P_mg_kgchar_unfilt_Mean ~ 1 + Solid_P_g_kg + Mono + Di + XANES_Pi_Ca + XANES_Pi_Mg + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp
## Leachate_P_mg_kgchar_unfilt_Mean ~ 1 + Solid_P_g_kg + Mono + 
##     Di + XANES_Pi_Ca + XANES_Pi_Mg + Burn_Severity + Land_Coverage_Category + 
##     Burn_Duration + Char_Max_Temp
Solid_P_g_kg ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp
## Solid_P_g_kg ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + 
##     Char_Max_Temp
Mono ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp
## Mono ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + 
##     Char_Max_Temp
Di ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp
## Di ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + 
##     Char_Max_Temp
XANES_Pi_Ca ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp
## XANES_Pi_Ca ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + 
##     Char_Max_Temp
XANES_Pi_Mg ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + Char_Max_Temp
## XANES_Pi_Mg ~ 1 + Burn_Severity + Land_Coverage_Category + Burn_Duration + 
##     Char_Max_Temp
model<-'
Solid_P_g_kg ~ 1 + Burn_Severity
Mono ~ 1 + Burn_Severity
Di ~ 1 + Burn_Severity
XANES_Pi_Ca ~ 1 + Burn_Severity
XANES_Pi_Mg ~ 1 + Burn_Severity 
Burn_Severity ~ 1 + Land_Coverage_Category + Burn_Duration + Char_Max_Temp
Char_Max_Temp ~ Land_Coverage_Category + Burn_Duration
'

fit <- sem(model, data = data_rda)
## Warning in lav_data_full(data = data, group = group, cluster = cluster, :
## lavaan WARNING: some observed variances are (at least) a factor 1000 times
## larger than others; use varTable(fit) to investigate
summary(fit, fit.measures=TRUE)
## lavaan 0.6.16 ended normally after 181 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        34
## 
##   Number of observations                            18
## 
## Model Test User Model:
##                                                       
##   Test statistic                                63.543
##   Degrees of freedom                                15
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                               185.488
##   Degrees of freedom                                35
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.677
##   Tucker-Lewis Index (TLI)                       0.247
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)               -418.634
##   Loglikelihood unrestricted model (H1)             NA
##                                                       
##   Akaike (AIC)                                 905.268
##   Bayesian (BIC)                               935.541
##   Sample-size adjusted Bayesian (SABIC)        831.069
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.424
##   90 Percent confidence interval - lower         0.320
##   90 Percent confidence interval - upper         0.534
##   P-value H_0: RMSEA <= 0.050                    0.000
##   P-value H_0: RMSEA >= 0.080                    1.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.158
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate    Std.Err  z-value  P(>|z|)
##   Solid_P_g_kg ~                                        
##     Burn_Severity      -0.684    0.968   -0.707    0.480
##   Mono ~                                                
##     Burn_Severity       7.760    2.295    3.381    0.001
##   Di ~                                                  
##     Burn_Severity       2.949    0.715    4.126    0.000
##   XANES_Pi_Ca ~                                         
##     Burn_Severity     -11.707    2.779   -4.212    0.000
##   XANES_Pi_Mg ~                                         
##     Burn_Severity      -4.667    4.299   -1.086    0.278
##   Burn_Severity ~                                       
##     Lnd_Cvrg_Ctgry      0.042    0.412    0.103    0.918
##     Burn_Duration      -0.001    0.001   -1.077    0.282
##     Char_Max_Temp      -0.002    0.001   -2.095    0.036
##   Char_Max_Temp ~                                       
##     Lnd_Cvrg_Ctgry     29.882   99.784    0.299    0.765
##     Burn_Duration       0.501    0.127    3.941    0.000
## 
## Covariances:
##                    Estimate    Std.Err  z-value  P(>|z|)
##  .Solid_P_g_kg ~~                                       
##    .Mono              -21.931   10.877   -2.016    0.044
##    .Di                 -4.552    3.167   -1.437    0.151
##    .XANES_Pi_Ca        21.671   12.665    1.711    0.087
##    .XANES_Pi_Mg        54.259   22.020    2.464    0.014
##  .Mono ~~                                               
##    .Di                 23.678    9.005    2.629    0.009
##    .XANES_Pi_Ca       -92.587   35.092   -2.638    0.008
##    .XANES_Pi_Mg       -95.291   48.076   -1.982    0.047
##  .Di ~~                                                 
##    .XANES_Pi_Ca       -27.776   10.774   -2.578    0.010
##    .XANES_Pi_Mg       -25.395   14.526   -1.748    0.080
##  .XANES_Pi_Ca ~~                                        
##    .XANES_Pi_Mg       122.063   58.965    2.070    0.038
## 
## Intercepts:
##                    Estimate    Std.Err  z-value  P(>|z|)
##    .Solid_P_g_kg        6.107    2.509    2.434    0.015
##    .Mono               -8.444    5.951   -1.419    0.156
##    .Di                 -4.745    1.853   -2.561    0.010
##    .XANES_Pi_Ca        64.778    7.206    8.990    0.000
##    .XANES_Pi_Mg        31.576   11.146    2.833    0.005
##    .Burn_Severity       3.592    0.698    5.146    0.000
##    .Char_Max_Temp     254.362  158.596    1.604    0.109
## 
## Variances:
##                    Estimate    Std.Err  z-value  P(>|z|)
##    .Solid_P_g_kg       17.123    5.708    3.000    0.003
##    .Mono               96.287   32.096    3.000    0.003
##    .Di                  9.336    3.112    3.000    0.003
##    .XANES_Pi_Ca       141.180   47.060    3.000    0.003
##    .XANES_Pi_Mg       337.761  112.587    3.000    0.003
##    .Burn_Severity       0.502    0.167    3.000    0.003
##    .Char_Max_Temp   29645.894 9881.965    3.000    0.003
modindices(fit, sort=TRUE)
##                        lhs op                    rhs     mi     epc sepc.lv
## 94           Char_Max_Temp  ~            XANES_Pi_Ca 10.746  11.939  11.939
## 93           Char_Max_Temp  ~                     Di  9.911 -44.588 -44.588
## 92           Char_Max_Temp  ~                   Mono  9.856 -13.846 -13.846
## 87           Burn_Severity  ~                   Mono  8.165  -0.068  -0.068
## 97  Land_Coverage_Category  ~           Solid_P_g_kg  7.567   0.064   0.064
## 109          Burn_Duration  ~            XANES_Pi_Mg  7.452  13.558  13.558
## 89           Burn_Severity  ~            XANES_Pi_Ca  7.164   0.053   0.053
## 101 Land_Coverage_Category  ~            XANES_Pi_Mg  5.624   0.012   0.012
## 85             XANES_Pi_Mg  ~          Burn_Duration  5.609   0.024   0.024
## 90           Burn_Severity  ~            XANES_Pi_Mg  5.516   0.030   0.030
## 56            Solid_P_g_kg  ~ Land_Coverage_Category  5.464   3.443   3.443
## 88           Burn_Severity  ~                     Di  5.337  -0.178  -0.178
## 40            Solid_P_g_kg ~~          Burn_Severity  3.020  -1.106  -1.106
## 48             XANES_Pi_Mg ~~          Burn_Severity  2.986   4.796   4.796
## 57            Solid_P_g_kg  ~          Burn_Duration  2.968  -0.004  -0.004
## 105          Burn_Duration  ~           Solid_P_g_kg  2.585  35.466  35.466
## 63                    Mono  ~ Land_Coverage_Category  2.413   4.119   4.119
## 42                    Mono ~~          Burn_Severity  2.370  -1.764  -1.764
## 62                    Mono  ~          Char_Max_Temp  1.955  -0.009  -0.009
## 55            Solid_P_g_kg  ~          Char_Max_Temp  1.681  -0.005  -0.005
## 100 Land_Coverage_Category  ~            XANES_Pi_Ca  1.595   0.010   0.010
## 64                    Mono  ~          Burn_Duration  1.315  -0.005  -0.005
## 108          Burn_Duration  ~            XANES_Pi_Ca  1.223   8.496   8.496
## 106          Burn_Duration  ~                   Mono  1.149  -9.970  -9.970
## 83             XANES_Pi_Mg  ~          Char_Max_Temp  1.010   0.016   0.016
## 77             XANES_Pi_Ca  ~ Land_Coverage_Category  0.897   3.312   3.312
## 95           Char_Max_Temp  ~            XANES_Pi_Mg  0.840   2.158   2.158
## 91           Char_Max_Temp  ~           Solid_P_g_kg  0.824   9.496   9.496
## 86           Burn_Severity  ~           Solid_P_g_kg  0.754   0.049   0.049
## 71                      Di  ~          Burn_Duration  0.673   0.001   0.001
## 76             XANES_Pi_Ca  ~          Char_Max_Temp  0.537   0.006   0.006
## 98  Land_Coverage_Category  ~                   Mono  0.500  -0.007  -0.007
## 99  Land_Coverage_Category  ~                     Di  0.404  -0.020  -0.020
## 78             XANES_Pi_Ca  ~          Burn_Duration  0.341  -0.003  -0.003
## 45                      Di ~~          Char_Max_Temp  0.340 -43.373 -43.373
## 107          Burn_Duration  ~                     Di  0.297 -16.287 -16.287
## 44                      Di ~~          Burn_Severity  0.174   0.169   0.169
## 84             XANES_Pi_Mg  ~ Land_Coverage_Category  0.109  -2.117  -2.117
## 70                      Di  ~ Land_Coverage_Category  0.109  -0.308  -0.308
## 46             XANES_Pi_Ca ~~          Burn_Severity  0.059   0.368   0.368
## 69                      Di  ~          Char_Max_Temp  0.007   0.000   0.000
##     sepc.all sepc.nox
## 94     0.824    0.824
## 93    -0.783   -0.783
## 92    -0.716   -0.716
## 87    -0.852   -0.852
## 97     0.597    0.597
## 109    0.732    0.732
## 89     0.880    0.880
## 101    0.524    0.524
## 85     0.443    0.001
## 90     0.566    0.566
## 56     0.368    0.821
## 88    -0.752   -0.752
## 40    -0.377   -0.377
## 48     0.368    0.368
## 57    -0.334   -0.001
## 105    0.423    0.423
## 63     0.147    0.328
## 42    -0.254   -0.254
## 62    -0.175   -0.175
## 55    -0.270   -0.270
## 100    0.381    0.381
## 64    -0.134    0.000
## 108    0.405    0.405
## 106   -0.356   -0.356
## 83     0.202    0.202
## 77     0.089    0.198
## 95     0.169    0.169
## 91     0.164    0.164
## 86     0.205    0.205
## 71     0.099    0.000
## 76     0.091    0.091
## 98    -0.194   -0.194
## 99    -0.190   -0.190
## 78    -0.067    0.000
## 45    -0.082   -0.082
## 107   -0.197   -0.197
## 44     0.078    0.078
## 84    -0.050   -0.112
## 70    -0.032   -0.072
## 46     0.044    0.044
## 69     0.011    0.011
model<-'
Leachate_P_mg_kgchar_unfilt_Mean ~ 1 + Solid_P_g_kg + Mono + Di + XANES_Pi_Ca + XANES_Pi_Mg + Burn_Severity
Solid_P_g_kg ~ 1 + Burn_Severity
Mono ~ 1 + Burn_Severity
Di ~ 1 + Burn_Severity
XANES_Pi_Ca ~ 1 + Burn_Severity
XANES_Pi_Mg ~ 1 + Burn_Severity 
'

fit <- sem(model, data = data_rda)
## Warning in lav_data_full(data = data, group = group, cluster = cluster, :
## lavaan WARNING: some observed variances are (at least) a factor 1000 times
## larger than others; use varTable(fit) to investigate
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some observed variances are larger than 1000000
##   lavaan NOTE: use varTable(fit) to investigate
summary(fit, fit.measures=TRUE)
## lavaan 0.6.16 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        23
## 
##   Number of observations                            18
## 
## Model Test User Model:
##                                                       
##   Test statistic                                60.715
##   Degrees of freedom                                10
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                               154.464
##   Degrees of freedom                                21
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.620
##   Tucker-Lewis Index (TLI)                       0.202
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)               -464.191
##   Loglikelihood unrestricted model (H1)             NA
##                                                       
##   Akaike (AIC)                                 974.382
##   Bayesian (BIC)                               994.861
##   Sample-size adjusted Bayesian (SABIC)        924.189
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.531
##   90 Percent confidence interval - lower         0.407
##   90 Percent confidence interval - upper         0.663
##   P-value H_0: RMSEA <= 0.050                    0.000
##   P-value H_0: RMSEA >= 0.080                    1.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.267
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                                      Estimate     Std.Err    z-value  P(>|z|)
##   Leachate_P_mg_kgchar_unfilt_Mean ~                                         
##     Solid_P_g_kg                         748.443     66.843   11.197    0.000
##     Mono                                  18.575     28.188    0.659    0.510
##     Di                                   -98.463     90.523   -1.088    0.277
##     XANES_Pi_Ca                          -30.447     23.279   -1.308    0.191
##     XANES_Pi_Mg                          -13.373     15.050   -0.889    0.374
##     Burn_Severity                         11.290    525.112    0.021    0.983
##   Solid_P_g_kg ~                                                             
##     Burn_Severity                         -0.684      0.968   -0.707    0.480
##   Mono ~                                                                     
##     Burn_Severity                          7.760      2.295    3.381    0.001
##   Di ~                                                                       
##     Burn_Severity                          2.949      0.715    4.126    0.000
##   XANES_Pi_Ca ~                                                              
##     Burn_Severity                        -11.707      2.779   -4.212    0.000
##   XANES_Pi_Mg ~                                                              
##     Burn_Severity                         -4.667      4.299   -1.086    0.278
## 
## Intercepts:
##                    Estimate     Std.Err    z-value  P(>|z|)
##    .Lcht_P_mg_k__M    -213.388   1847.683   -0.115    0.908
##    .Solid_P_g_kg         6.107      2.509    2.434    0.015
##    .Mono                -8.444      5.951   -1.419    0.156
##    .Di                  -4.745      1.853   -2.561    0.010
##    .XANES_Pi_Ca         64.778      7.206    8.990    0.000
##    .XANES_Pi_Mg         31.576     11.146    2.833    0.005
## 
## Variances:
##                    Estimate     Std.Err    z-value  P(>|z|)
##    .Lcht_P_mg_k__M 1377085.799 459028.600    3.000    0.003
##    .Solid_P_g_kg        17.123      5.708    3.000    0.003
##    .Mono                96.287     32.096    3.000    0.003
##    .Di                   9.336      3.112    3.000    0.003
##    .XANES_Pi_Ca        141.180     47.060    3.000    0.003
##    .XANES_Pi_Mg        337.761    112.587    3.000    0.003
modindices(fit, sort=TRUE)
##              lhs op                              rhs     mi     epc sepc.lv
## 36          Mono ~~                      XANES_Pi_Ca 11.351 -92.587 -92.587
## 58   XANES_Pi_Ca  ~                             Mono 11.351  -0.962  -0.962
## 49          Mono  ~                      XANES_Pi_Ca 11.351  -0.656  -0.656
## 35          Mono ~~                               Di 11.226  23.678  23.678
## 53            Di  ~                             Mono 11.226   0.246   0.246
## 48          Mono  ~                               Di 11.226   2.536   2.536
## 38            Di ~~                      XANES_Pi_Ca 10.536 -27.776 -27.776
## 59   XANES_Pi_Ca  ~                               Di 10.536  -2.975  -2.975
## 54            Di  ~                      XANES_Pi_Ca 10.536  -0.197  -0.197
## 45  Solid_P_g_kg  ~                      XANES_Pi_Mg  9.163   0.161   0.161
## 34  Solid_P_g_kg ~~                      XANES_Pi_Mg  9.163  54.259  54.259
## 62   XANES_Pi_Mg  ~                     Solid_P_g_kg  9.163   3.169   3.169
## 61   XANES_Pi_Mg  ~ Leachate_P_mg_kgchar_unfilt_Mean  7.661   0.004   0.004
## 41  Solid_P_g_kg  ~ Leachate_P_mg_kgchar_unfilt_Mean  6.033  -0.004  -0.004
## 65   XANES_Pi_Mg  ~                      XANES_Pi_Ca  5.624   0.865   0.865
## 60   XANES_Pi_Ca  ~                      XANES_Pi_Mg  5.624   0.361   0.361
## 40   XANES_Pi_Ca ~~                      XANES_Pi_Mg  5.624 122.063 122.063
## 47          Mono  ~                     Solid_P_g_kg  5.251  -1.281  -1.281
## 31  Solid_P_g_kg ~~                             Mono  5.251 -21.931 -21.931
## 42  Solid_P_g_kg  ~                             Mono  5.251  -0.228  -0.228
## 63   XANES_Pi_Mg  ~                             Mono  5.026  -0.990  -0.990
## 37          Mono ~~                      XANES_Pi_Mg  5.026 -95.291 -95.291
## 50          Mono  ~                      XANES_Pi_Mg  5.026  -0.282  -0.282
## 46          Mono  ~ Leachate_P_mg_kgchar_unfilt_Mean  4.065  -0.001  -0.001
## 64   XANES_Pi_Mg  ~                               Di  3.681  -2.720  -2.720
## 39            Di ~~                      XANES_Pi_Mg  3.681 -25.395 -25.395
## 55            Di  ~                      XANES_Pi_Mg  3.681  -0.075  -0.075
## 33  Solid_P_g_kg ~~                      XANES_Pi_Ca  3.497  21.671  21.671
## 44  Solid_P_g_kg  ~                      XANES_Pi_Ca  3.497   0.153   0.153
## 57   XANES_Pi_Ca  ~                     Solid_P_g_kg  3.497   1.266   1.266
## 56   XANES_Pi_Ca  ~ Leachate_P_mg_kgchar_unfilt_Mean  3.175   0.002   0.002
## 32  Solid_P_g_kg ~~                               Di  2.333  -4.552  -4.552
## 43  Solid_P_g_kg  ~                               Di  2.333  -0.488  -0.488
## 52            Di  ~                     Solid_P_g_kg  2.333  -0.266  -0.266
## 51            Di  ~ Leachate_P_mg_kgchar_unfilt_Mean  0.623   0.000   0.000
## 66 Burn_Severity  ~ Leachate_P_mg_kgchar_unfilt_Mean  0.000   0.000   0.000
## 67 Burn_Severity  ~                     Solid_P_g_kg  0.000   0.000   0.000
## 71 Burn_Severity  ~                      XANES_Pi_Mg  0.000   0.000   0.000
## 24 Burn_Severity ~~                    Burn_Severity  0.000   0.000   0.000
## 70 Burn_Severity  ~                      XANES_Pi_Ca  0.000   0.000   0.000
## 68 Burn_Severity  ~                             Mono  0.000   0.000   0.000
## 69 Burn_Severity  ~                               Di  0.000   0.000   0.000
##    sepc.all sepc.nox
## 36   -0.794   -0.794
## 58   -0.721   -0.721
## 49   -0.875   -0.875
## 35    0.790    0.790
## 53    0.724    0.724
## 48    0.862    0.862
## 38   -0.765   -0.765
## 59   -0.757   -0.757
## 54   -0.773   -0.773
## 45    0.726    0.726
## 34    0.713    0.713
## 62    0.701    0.701
## 61    0.678    0.678
## 41   -3.426   -3.426
## 65    0.763    0.763
## 60    0.409    0.409
## 40    0.559    0.559
## 47   -0.428   -0.428
## 31   -0.540   -0.540
## 42   -0.681   -0.681
## 63   -0.655   -0.655
## 37   -0.528   -0.528
## 50   -0.427   -0.427
## 46   -0.398   -0.398
## 64   -0.611   -0.611
## 39   -0.452   -0.452
## 55   -0.335   -0.335
## 33    0.441    0.441
## 44    0.613    0.613
## 57    0.317    0.317
## 56    0.321    0.321
## 32   -0.360   -0.360
## 43   -0.495   -0.495
## 52   -0.262   -0.262
## 51   -0.143   -0.143
## 66    0.000    0.000
## 67    0.000    0.000
## 71    0.000    0.000
## 24    0.000    0.000
## 70    0.000    0.000
## 68    0.000    0.000
## 69    0.000    0.000
#reference: https://kevintshoemaker.github.io/NRES-746/SEM.RMarkdown.html
#making a structural regression
model<-'
chem =~ Solid_P_g_kg + Mono + Di + XANES_Pi_Ca + XANES_Pi_Mg 
mobilization =~ Leachate_P_mg_kgchar_filt0.7_Mean + MBD_P_mg_kgchar_Mean
mobilization ~ chem + Burn_Severity
'

fit<-sem(model, data = data_rda)
## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov
## variances are negative
summary(fit, fit.measures=TRUE, standardized=TRUE)
## lavaan 0.6.16 ended normally after 335 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        16
## 
##   Number of observations                            18
## 
## Model Test User Model:
##                                                       
##   Test statistic                                41.992
##   Degrees of freedom                                19
##   P-value (Chi-square)                           0.002
## 
## Model Test Baseline Model:
## 
##   Test statistic                               156.638
##   Degrees of freedom                                28
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.821
##   Tucker-Lewis Index (TLI)                       0.737
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)               -488.780
##   Loglikelihood unrestricted model (H1)             NA
##                                                       
##   Akaike (AIC)                                1009.559
##   Bayesian (BIC)                              1023.805
##   Sample-size adjusted Bayesian (SABIC)        974.642
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.259
##   90 Percent confidence interval - lower         0.153
##   90 Percent confidence interval - upper         0.366
##   P-value H_0: RMSEA <= 0.050                    0.003
##   P-value H_0: RMSEA >= 0.080                    0.994
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.226
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate   Std.Err   z-value  P(>|z|)   Std.lv   Std.all
##   chem =~                                                                  
##     Solid_P_g_kg       1.000                                 2.066    0.493
##     Mono              -5.738     2.494   -2.301    0.021   -11.856   -0.945
##     Di                -1.911     0.836   -2.286    0.022    -3.948   -0.926
##     XANES_Pi_Ca        7.601     3.311    2.295    0.022    15.706    0.938
##     XANES_Pi_Mg        5.381     2.910    1.849    0.064    11.120    0.586
##   mobilization =~                                                          
##     Lcht_P___0.7_M     1.000                                98.360    0.829
##     MBD_P_mg_kgc_M     1.014     2.044    0.496    0.620    99.733    1.158
## 
## Regressions:
##                    Estimate   Std.Err   z-value  P(>|z|)   Std.lv   Std.all
##   mobilization ~                                                           
##     chem               0.662     8.050    0.082    0.934     0.014    0.014
##     Burn_Severity     -4.959    27.075   -0.183    0.855    -0.050   -0.051
## 
## Variances:
##                    Estimate   Std.Err   z-value  P(>|z|)   Std.lv   Std.all
##    .Solid_P_g_kg      13.328     4.507    2.957    0.003    13.328    0.757
##    .Mono              16.870     9.445    1.786    0.074    16.870    0.107
##    .Di                 2.581     1.212    2.129    0.033     2.581    0.142
##    .XANES_Pi_Ca       33.668    17.467    1.928    0.054    33.668    0.120
##    .XANES_Pi_Mg      236.213    80.621    2.930    0.003   236.213    0.656
##    .Lcht_P___0.7_M  4417.761 19515.513    0.226    0.821  4417.761    0.313
##    .MBD_P_mg_kgc_M -2529.317 20024.632   -0.126    0.899 -2529.317   -0.341
##     chem               4.270     3.905    1.093    0.274     1.000    1.000
##    .mobilization    9647.787 19723.108    0.489    0.625     0.997    0.997
# Define the SEM model
model <- '
  # Measurement Model
    Chem =~ Solid_P_g_kg + Mono + Di + XANES_Pi_Ca + XANES_Pi_Mg
    Mobilization =~ Leachate_P_mg_kgchar_unfilt_Mean + MBD_P_mg_kgchar_Mean

  # Structural Model
    Mobilization ~ c*Burn_Severity + d*Chem + e*Burn_Severity*Chem

  # Residual Variances
    Land_Coverage_Category ~~ Land_Coverage_Category
    Burn_Duration ~~ Burn_Duration
    Char_Max_Temp ~~ Char_Max_Temp
    Solid_P_g_kg ~~ Solid_P_g_kg
    Mono ~~ Mono
    Di ~~ Di
    XANES_Pi_Ca ~~ XANES_Pi_Ca
    XANES_Pi_Mg ~~ XANES_Pi_Mg
    Leachate_P_mg_kgchar_unfilt_Mean ~~ Leachate_P_mg_kgchar_unfilt_Mean
    MBD_P_mg_kgchar_Mean ~~ MBD_P_mg_kgchar_Mean
'

# Fit the SEM model
fit <- sem(model, data = data_rda)
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some observed variances are (at least) a factor 1000 times larger than others; use varTable(fit) to investigate

## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: some observed variances are larger than 1000000
##   lavaan NOTE: use varTable(fit) to investigate
## Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
##     Could not compute standard errors! The information matrix could
##     not be inverted. This may be a symptom that the model is not
##     identified.
# Summarize the results
summary(fit, fit.measures=TRUE, standardized=TRUE)
## lavaan 0.6.16 ended normally after 119 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        19
## 
##   Number of observations                            18
## 
## Model Test User Model:
##                                                       
##   Test statistic                               175.445
##   Degrees of freedom                                46
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                               273.519
##   Degrees of freedom                                55
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.408
##   Tucker-Lewis Index (TLI)                       0.292
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)               -834.409
##   Loglikelihood unrestricted model (H1)             NA
##                                                       
##   Akaike (AIC)                                1706.819
##   Bayesian (BIC)                              1723.736
##   Sample-size adjusted Bayesian (SABIC)       1665.355
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.395
##   90 Percent confidence interval - lower         0.334
##   90 Percent confidence interval - upper         0.458
##   P-value H_0: RMSEA <= 0.050                    0.000
##   P-value H_0: RMSEA >= 0.080                    1.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.367
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate     Std.Err  z-value  P(>|z|)   Std.lv     Std.all
##   Chem =~                                                                     
##     Solid_P_g_kg         1.000                                  2.188    0.522
##     Mono                -5.444       NA                       -11.914   -0.950
##     Di                  -1.793       NA                        -3.924   -0.921
##     XANES_Pi_Ca          7.126       NA                        15.596    0.931
##     XANES_Pi_Mg          5.188       NA                        11.353    0.598
##   Mobilization =~                                                             
##     Lcht_P_mg_k__M       1.000                               2123.611    0.619
##     MBD_P_mg_kgc_M      -0.005       NA                       -10.906   -0.126
## 
## Regressions:
##                    Estimate     Std.Err  z-value  P(>|z|)   Std.lv     Std.all
##   Mobilization ~                                                              
##     Burn_Svrty (c)    1005.436       NA                         0.473    0.477
##     Chem       (d)     852.823       NA                         0.879    0.879
## 
## Variances:
##                    Estimate     Std.Err  z-value  P(>|z|)   Std.lv     Std.all
##     Lnd_Cvrg_Ctgry       0.201       NA                         0.201    1.000
##     Burn_Duration   123577.396       NA                    123577.396    1.000
##     Char_Max_Temp    58868.714       NA                     58868.714    1.000
##    .Solid_P_g_kg        12.809       NA                        12.809    0.728
##    .Mono                15.484       NA                        15.484    0.098
##    .Di                   2.767       NA                         2.767    0.152
##    .XANES_Pi_Ca         37.119       NA                        37.119    0.132
##    .XANES_Pi_Mg        230.993       NA                       230.993    0.642
##    .Lcht_P_mg_k__M 7276330.946       NA                   7276330.946    0.617
##    .MBD_P_mg_kgc_M    7346.177       NA                      7346.177    0.984
##     Chem                 4.789       NA                         1.000    1.000
##    .Mobilization         0.678       NA                         0.000    0.000

#haven’t updated below here yet (NMR and XANES code) fix

#spectra stacked #spectra <- read.csv(“~/Library/CloudStorage/OneDrive-SharedLibraries-PNNL/RC-3, River Corridor SFA - Presentations/ESA_2023/Morgan/spectra.csv”)

#filt<-spectra[ -c(6:9)] #filt\(P_BSLE_0007_filt07_sample<-filt\)P_BSLE_0007_filt07_sample+7 #filt\(P_BSLE_0002_filt07_sample<-filt\)P_BSLE_0002_filt07_sample+14 #filt\(P_BSLE_0050_filt07_sample<-filt\)P_BSLE_0050_filt07_sample+21 #filt<-melt(filt, id=c(“energy”))

#ggplot(filt, aes(x=energy, y=value))+geom_line(aes(color = variable), show.legend = FALSE)+theme_classic()+ylab(“Noramlized Absorption (arb. units)”)+ xlab(“Energy (eV)”)+theme(text = element_text(size = 15))+xlim(2140,2200)+scale_color_manual(values=c(‘Black’,‘Black’, ‘Black’, ‘Black’))

#plot spectra of samples #BSLE_solid_sample_spectra_normalizedforfinalLCF_12.1.22 <- read.csv(“~/PNNL/RC-3, River Corridor SFA - Burn Severity Experiments/Burn Severity Experiment/Data/P_all data streams/P XAS/Processed_final/BSLE_Solids_Final_Fits/BSLE_solid_sample_spectra_normalizedforfinalLCF_12.1.22.csv”)

#DFL<-BSLE_solid_sample_spectra_normalizedforfinalLCF_12.1.22[ c(25, 20, 19, 13, 1)] #DFL<-melt(DFL, id=c(“energy”))

#ggplot(DFL, aes(x=energy, y=value))+geom_line()+theme_classic()+ylab(“Noramlized Absorption (arb. units)”)+ xlab(“Energy (eV)”)+facet_grid(variable~.)+theme(text = element_text(size = 15))+xlim(2140,2200)

#try to stack them differently #DFL<-BSLE_solid_sample_spectra_normalizedforfinalLCF_12.1.22[ c(25, 20, 19, 13, 1)] #DFL\(P_BSLE_0073_solidchar_036_fl_4s_avg<-DFL\)P_BSLE_0073_solidchar_036_fl_4s_avg+9 #DFL\(P_BSLE_0052_solidchar_035_fl_6sc_avg<-DFL\)P_BSLE_0052_solidchar_035_fl_6sc_avg+6 #DFL\(P_BSLE_0051_solidchar_034_fl_6s_avg<-DFL\)P_BSLE_0051_solidchar_034_fl_6s_avg+3 #DFL<-melt(DFL, id=c(“energy”))

#DFL.sample.xas.spectra.gradient<-ggplot(DFL, aes(x=energy, y=value))+geom_line(aes(color = variable), show.legend = FALSE)+theme_classic()+ylab(“Noramlized Absorption (arb. units)”)+ xlab(“Energy (eV)”)+theme(text = element_text(size = 15))+xlim(2140,2200)+scale_color_manual(values=c(‘Black’,‘Black’, ‘Black’, ‘Black’))